domain: info math spi

text:
Abla Kammoun

Axel Müller

Emil Björnson

Mérouane Debbah

LOW-COMPLEXITY LINEAR PRECODING FOR MULTI-CELL MASSIVE MIMO SYSTEMS

Keywords: Massive MIMO, linear precoding, low complexity, multi-cell systems, random matrix theory

HAL is

INTRODUCTION

One of the major challenges in multi-cell systems is dealing with inter-cell and intra-cell interference. Currently, the general trend is to deploy multiple antennas at the base stations (BSs), thereby enabling multi-user MIMO with flexible spatial interference mitigation [START_REF] Gesbert | Multi-cell MIMO cooperative networks: A new look at interference[END_REF]. User separation in the downlink is then performed using linear precoding. Unfortunately, the use of optimal precoding designs is far from being incorporated in current wireless standards such as LTE-Advanced [START_REF] Holma | LTE advanced: 3GPP solution for IMT-Advanced[END_REF]. This can be attributed to the fact that very accurate instantaneous channel state information (CSI) is required, which can be cumbersome to achieve in practice [START_REF] Björnson | Optimal resource allocation in coordinated multi-cell systems[END_REF]. One can imagine that this becomes worse as the number of BS antennas, M , and the number of users, K, increase. Interestingly, this is not the case in reality. As M → ∞ for a This research has been supported by the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for Complex Network Engineering) and an International Postdoc Grant from The Swedish Research Council. fixed K, simple linear precoding, like maximum ratio transmission (MRT), is asymptotically optimal [START_REF] Rusek | Scaling up MIMO: Opportunities and Challenges with Very Large Arrays[END_REF] and robust to CSI imperfections [START_REF] Marzetta | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF]. Nevertheless, MRT is not very appealing at practical values of M and K, since it does not actively suppress residual inter-user interference [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF]. In fact, a precoding design based on both M and K growing large, with a fixed ratio, yields better massive MIMO performance [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF]. In this regime, RZF precoding is near-optimal from a throughput perspective. However, it requires calculation of the inverse of the Gram matrix of the joint channel of all users, which is a non-trivial operation with a complexity scaling of M K 2 . Fortunately, the system performance is predictable in the large-(M, K) regime, where advanced tools from random matrix theory provide deterministic approximations of the achievable rates [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF]. In light of these results, we proposed in [START_REF] Müller | Linear precoding based on truncated polynomial expansion-Part I: Large-scale single-cell systems[END_REF] to solve the precoding complexity issues in single-cell systems using a new family of precoding schemes called TPE precoding. This family was obtained by approximating the matrix inverse in RZF precoding by a matrix polynomial. A similar approach was independently proposed in [START_REF] Zarei | Lowcomplexity linear precoding for downlink large-scale MIMO systems[END_REF].

In this paper, we extend the TPE precoding from [START_REF] Müller | Linear precoding based on truncated polynomial expansion-Part I: Large-scale single-cell systems[END_REF] to the scenario of multi-cell massive MIMO systems. A special focus is placed on the analysis of realistic characteristics, including user-specific channel covariance matrices, imperfect CSI and pilot contamination. The proofs of our results can be found in the extended version of this paper; see [START_REF] Kammoun | Linear Precoding Based on Truncated Polynomial Expansion-Part II: Large-Scale Multi-Cell Systems[END_REF].

SYSTEM MODEL

We consider the downlink of a multi-cell system composed of L > 1 cells. Each cell consists of a M -antenna BS serving K single-antenna users. We assume a time-division duplex (TDD) protocol where the BS acquires instantaneous CSI in the uplink and uses it for the downlink transmission, exploiting channel reciprocity. The TDD protocols are synchronized across cells, so that pilot signaling and data transmission take place simultaneously. The received complex baseband signal at the mth user terminal (UT) in the jth cell is

y j,m = L ℓ=1 h H ℓ,j,m x ℓ + n j,m (1)
where x ℓ ∈ C M ×1 is the transmit signal from the ℓth BS and h ℓ,j,m ∈ C M ×1 is the channel vector from that BS to the mth UT in the jth cell, and n j,m ∼ CN (0, σ 2 ) is the additive circularly-symmetric complex Gaussian noise with variance σ 2 . The channel vectors are modeled as Rayleigh fading with

h ℓ,j,m ∼ CN 0, R ℓ,j,m (2)
where the family of covariance matrices (R ℓ,j,m ) L,L,K ℓ=1,j=1,m=1

satisfy the following conditions:

• Bounded norm: lim sup M R ℓ,j,m 2 < +∞, ∀ℓ, j, m;

• Trace scaling: lim inf M 1 M tr (R ℓ,j,m ) > 0, ∀ℓ, j, m; • Finite dimensional matrix space for R ℓ,j,m : It exists a finite integer S > 0 and a linear independent family of matrices F 1 , . . . , F S such that R ℓ,j,m = S k=1 α ℓ,j,m,k F k . Note that these conditions are less restrictive than the one used in [START_REF] Hoydis | Asymptotic moments for interference mitigation in correlated fading channels[END_REF], where R ℓ,j,m was assumed to belong to a finite set of matrices. It is also in agreement with several physical channel models presented in the literature; for example, the one-ring model with user groups from [START_REF] Adhikary | Joint spatial division and multiplexing-the large-scale array regime[END_REF]. This channel model considers a finite number G of groups which share approximately the same location and thus the same covariance matrix. Let θ ℓ,j,g and ∆ ℓ,j,g be, respectively, the azimuth angle and the azimuth angular spread between the cell ℓ and the users in group g of cell j. Moreover, let d be the distance between two adjacent antennas (see Fig. 1 in [START_REF] Adhikary | Joint spatial division and multiplexing-the large-scale array regime[END_REF]). Then, the (u, v)th entry of the covariance matrix R ℓ,j,m for users in group g is

[R ℓ,j,m ] u,v = 1 2∆ ℓ,j,g ∆ ℓ,j,g +θ ℓ,j,g -∆ ℓ,j,g +θ ℓ,j,g e d(u-v) sin α dα. (3)
We also assume that all BSs use Gaussian codebooks and linear precoding, such that the jth cell transmits the signal

x j = K m=1 g j,m s j,m = G j s j (4)
where G j = [g j,1 , . . . , g j,K ] ∈ C M ×K is the precoding matrix and s j = [s j,1 , . . . , s j,K ] ∼ CN (0, I K ) is the vector containing the data symbols for UTs in the jth cell. The transmission at BS j is subject to a transmit power constraint

1 K tr G j G H j = P j (5)
where P j is the average transmit power per user in the jth cell. The received signal (1) can be thus expressed as

y j,m = L ℓ=1 K k=1 h H ℓ,j,m g ℓ,k s ℓ,k + n j,m . (6)
A well known feature of large-scale MIMO systems is channel hardening, which implies that the effective useful channel h H j,j,m g j,m converges to its average value as M → ∞. We decompose the received signal as

yj,m = E h H j,j,m gj,m sj,m + sj,m h H j,j,m gj,m -E h H j,j,m gj,m + (ℓ,k) =(j,m) h ℓ,j,m g ℓ,k s ℓ,k + nj,m
and assume, similar to [START_REF] Marzetta | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF][START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF], that the receiver only knows the average channel gain E[h H j,j,m g j,m ], the average sum interference power ℓ,k E[|h H ℓ,j,m g ℓ,k | 2 ] and the variance of the noise σ 2 . By treating interference as worst-case Gaussian noise, the ergodic achievable rate of UT m in cell j is

r j,m = log 2 (1 + γ j,m ) (7)
where

γ j,m = |E [h j,j,m g j,m ]| 2 σ 2 + ℓ,k E |h H ℓ,j,m g ℓ,k | 2 -E h H j,j,m g j,m
2 .

(8)

Model of Imperfect Channel State information

Based on a TDD protocol, the channel is estimated in the uplink. In each cell, the UTs transmit mutually orthogonal pilot sequences, thereby allowing the BS to acquire CSI. Since the same set of orthogonal sequences is reused in each cell, the channel estimation is corrupted by inter-cell interference; this is called pilot contamination [START_REF] Marzetta | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF]. To estimate the channel corresponding to UT k in cell j, each BS correlates the received signal with the pilot sequence of that user. This results in the processed received signal

y tr j,k = h j,j,k + ℓ =j h j,ℓ,k + 1 √ ρ tr n tr j,k (9)
where n tr j,k ∼ CN (0, I M ) and ρ tr > 0 is the effective pilot SNR [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF]. The MMSE estimate h j,j,k of h j,j,k is

h j,j,k = R j,j,k S j,k y tr j,k = R j,j,k S j,k L ℓ=1 h j,ℓ,k + 1 √ ρ tr n tr j,k where S j,k = ( 1 ρtr + L ℓ=1 R ℓ,j,k ) -1 .
The estimated channel vectors at the jth BS to all UTs in its cell is denoted by

H j,j = h j,j,1 , . . . , h j,j,K ∈ C M ×K . ( 10
)
We also define Φ j,ℓ,k = R j,j,k S j,k R j,ℓ,k ∈ C M ×M and note that h j,j,k ∼ CN (0, Φ j,j,k ) is independent from the estimation error h j,j,kh j,j,k since the MMSE estimator is used.

MULTI-CELL LINEAR PRECODING

Regularized zero-forcing precoding

For multi-cell systems, the optimal linear precoding is unknown under imperfect CSI and requires extensive optimization procedures under perfect CSI [START_REF] Björnson | Optimal resource allocation in coordinated multi-cell systems[END_REF]. Therefore, only heuristic precoding schemes are feasible in multi-cell systems. RZF precoding is the state-of-the-art heuristic scheme in terms of system throughput [START_REF] Björnson | Optimal resource allocation in coordinated multi-cell systems[END_REF]. Using the notation of [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF], the RZF precoding matrix used by the BS in the jth cell is

G rzf j = √ Kβ j H j,j H H j,j + Kϕ j I M -1 H j,j (11)
where the scalar parameter β j > 0 is set to satisfy the power constraint in (5) and ϕ j is a positive regularizing parameter.

Prior works have considered the optimization of the parameter ϕ j in the single-cell case. This parameter provides a balance between maximization of the channel gain at each intended receiver (ϕ j is large) and the suppression of inter-user interference (when ϕ j is small). To the authors' knowledge, a closed-form optimization of the regularization parameter for multi-cell scenarios has thus far not been achieved. Hence, previous works have been restricted to the analysis of intuitive choices of the regularization parameter ϕ j . In this context, the work in [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF] considers massive MIMO performance of the RZF precoding in the large-(M, K) regime, where M and K tend to infinity such that

0 < lim inf K M ≤ lim sup K M < +∞. (12)

Truncated Polynomial Expansion Precoding

The matrix multiplications and inversion in [START_REF] Adhikary | Joint spatial division and multiplexing-the large-scale array regime[END_REF] gives RZF precoding an unfavorable complexity scaling of M K 2 [START_REF] Müller | Linear precoding based on truncated polynomial expansion-Part I: Large-scale single-cell systems[END_REF].

Building on the concept of TPE used in our work in [START_REF] Müller | Linear precoding based on truncated polynomial expansion-Part I: Large-scale single-cell systems[END_REF], we propose a new class of low-complexity linear precoding schemes also for the multi-cell case. The proposed precoding originates from the Cayley-Hamilton theorem which states that the inverse of a M × M matrix can be written as a weighted sum of its first M powers. A simplified precoding can then be obtained by considering only the first matrix powers. Assume that the jth BS employs a truncation order J j , then the proposed TPE precoding matrix is given by

G TPE j = Jj -1 n=0 w n,j H j,j H H j,j K n H j,j √ K
where {w n,j , j = 0, . . . , J j -1} are the J j scalar coefficients that are employed by the jth BS and the normalization by √ K controls the energy of the precoding matrix. Note that while the RZF precoding has only a single regularization parameter, the proposed TPE precoding scheme offers a larger set of J j design parameters. These coefficients define a parametrized class of precoding schemes ranging from MRT (J j = 1) to RZF precoding, which is achieved at J j = min (M, K) by selecting w n,j from the characteristic polynomial of

ASYMPTOTIC PERFORMANCE ANALYSIS

We provide in this section an asymptotic performance analysis of the proposed TPE precoding. In particular, we show that in the large (M, K)-regime, the SINR experienced by the mth UT served by the jth cell can be approximated by a deterministic term which depends only on the channel statistics. Before stating our main results, we cast the SINR expression (8) in a simpler form. Let w j = w 0,j , . . . , w Jj -1,j

T and let a j,m ∈ C Jj ×1 and B ℓ,j,m ∈ C Jj ×Jj be given by

[a j,m ] n = h H j,j,m √ K V n,j h H j,j,m √ K , n ∈ [0, J j -1] [B ℓ,j,m ] n,p = 1 K h H ℓ,j,m V n+p+1,ℓ h ℓ,j,m n, p ∈ [0, J ℓ -1]
where V n,j = Hj,j H H j,j K n . Then, the SINR experienced by the mth UT in the jth cell is

γ j,m = E w H j a j,m 2
σ 2 K + L ℓ=1 E [w H ℓ B ℓ,j,m w ℓ ] -E w H j a j,m 2 .
(13) As a j,m and B ℓ,j,m are of finite dimensions, it suffices to determine an asymptotic approximation of the expected value of each of their elements. For that, we introduce the functionals

X j,m (t) = 1 K h H j,j,m Σ(t, j) h j,j,m Z ℓ,j,m (t) = 1 K h H ℓ,j,m Σ(t, ℓ)h ℓ,j,m
where

Σ(t, j) = t Hj,j H H j,j K + I M -1 . It is easy to see that [a j,m ] n = (-1) n n! X (n) j,m (14)
[B ℓ,j,m ] n,p = (-1) n+p+1 (n + p + 1)! Z (n+p+1) ℓ,j,m . ( 15
)
where

X (n) j,m and Z (n+p+1) ℓ,j,m
represent the derivatives of X j,m (t) and Z ℓ,j,m (t) at t = 0.

Theorem 1 Let X j,m (t) and Z ℓ,j,m (t) be

X j,m (t) = δ j,m (t) 1 + tδ j,m (t) Z ℓ,j,m (t) = 1 K tr (R ℓ,j,m T ℓ (t)) - t 1 K tr Φ ℓ,j,m T ℓ (t) 2 1 + tδ ℓ,m (t)
where for each j = 1, . . . , L, m = 1, . . . , K, δ j,m (t) are the unique positive solutions to the following system of equations:

δ j,m (t) = 1 K tr   Φ j,j,m 1 K K k=1 tΦ j,j,k 1 + tδ j,k (t) + I M -1   and T j (t) = 1 K K k=1 tΦ j,j,k 1+tδ j,k (t) + I M -1
. In the large-(M, K) regime we have

E [X j,m (t)] -X j,m (t) -------→ M,K→+∞ 0 E [Z ℓ,j,m (t)] -Z ℓ,j,m (t) -------→ M,K→+∞ 0.

Corollary 2

The following holds in the large-(M, K) regime:

E X (n) j,m -X (n) j,m -------→ M,K→+∞ 0, (16) E Z (n) ℓ,j,m -Z (n) ℓ,j,m -------→ M,K→+∞ 0 (17)
where X

(n) j,m and Z

(n) ℓ,j,m are the derivatives of X(t) and Z ℓ,j,m (t), respectively, at t = 0.

The deterministic quantities X K tr(Φ j,j,m ) and Z ℓ,j,m = 1 K tr(R ℓ,j,m ). We can compute the deterministic sequences

X (n) j,m and Z (n) ℓ,j,m as X (n) j,m = - n k=1 n k kX (k-1) j,m δ (n-k) j,m + δ (n) j,m Z (n) ℓ,jab,m = 1 K tr R ℓ,j,m T (n) ℓ - n k=0 k n k δ (n-k) l,m Z (k-1) ℓ,j,m + n k=0 k n k δ (n-k) ℓ,m 1 K tr R ℓ,j,m T (k-1) ℓ where δ (p) ℓ,k and T (p) ℓ
can be computed using the iterative algorithm in [START_REF] Hoydis | Asymptotic moments for interference mitigation in correlated fading channels[END_REF]. Substituting the deterministic equivalent of Theorem 1 into (14) and (15), we get the following result.

Corollary 3 Let a j,m and B ℓ,j,m be given by

[a j,m ] n = (-1) n n! X (n) j,m , B ℓ,j,m n,p = (-1) n+p+1 (n + p + 1)! Z (n+p+1) ℓ,j,m
.

Then, a j,m and B ℓ,j,m converge as

max ℓ,j,m E B ℓ,j,m -B ℓ,j,m , E [a j,m -a j,m ] -------→ M,K→+∞ 0.
in the large-(M, K) regime and the SINRs converge as

γ j,m -γ j,m -------→ M,K→+∞ 0 ( 18
)
where γ j,m = w H j aj,ma H j,m wj L ℓ=1 w H ℓ B ℓ,j,m w ℓ -w H j aj,ma H j,m wj .

SYSTEM PERFORMANCE OPTIMIZATION

In the previous section, we derived deterministic equivalents of the SINR at each UT in the multi-cell system as a function of the polynomial coefficients {w ℓ,j , ℓ ∈ [1, L] , j ∈ [0, J ℓ -1]} of the TPE precoding. These coefficients can be selected to maximize any system performance metric. Furthermore, these coefficients need to be scaled to satisfy the transmit power constraints

1 K tr G ℓ,TPE G H ℓ,TPE = w T ℓ C ℓ w ℓ = P ℓ (19)
where C ℓ is a J ℓ × J ℓ matrix with elements given by

[C ℓ ] n,m = 1 K tr H ℓ,ℓ H ℓ,ℓ K n+m+1 , n, m ∈ [0, J ℓ -1] .
(20) We would like to pre-optimize the weights offline, thus the weights should not depend on the instantaneous value of the channel, but only its statistics. To this end, we substitute (19) by its asymptotic approximation

w T ℓ C ℓ w ℓ = P ℓ (21)
where

C ℓ n,m = (-1) n+m+1 (n+m+1)! 1 K tr T (n+m+1) ℓ
. In this paper, the performance metric is weighted maxmin fairness. In other words, we maximize the minimal value of

log 2 (1+γj,m) νj,m
for some user-specific weights ν j,m . If we would like to mimic the performance of RZF, then ν j,m should be the rate that user m in cell j achieves asymptotically with RZF precoding. Using deterministic equivalents, the corresponding optimization problem is

max w1,...,w L min j∈[1,L] m∈[1,K] log 2 1+ w H j aj,ma H j,m wj L ℓ=1 w H ℓ B ℓ,j,m w ℓ -w H j aj,ma H j,m wj ν j,m subject to w ℓ C ℓ w ℓ = P ℓ , ℓ ∈ [1, L] . ( 22
)
This problem is non-convex, but very similar to the multi-cast beamforming problems analyzed in [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF]. In particular, we can instead solve the following tractable relaxed convex problem:

subject to W ℓ 0, tr

C ℓ W ℓ = P ℓ , ℓ ∈ [1, L] a H j,m W j a j,m L ℓ=1 tr B ℓ,j,m W ℓ -a H j,m W j a j,m ≥ 2 νj.mξ -1, ∀j, m.

SIMULATION RESULTS

In this section, we consider the three-site sector in Fig. 1. Each of the three BSs is equipped with an horizontal linear array of M antennas. The radiation pattern of each antenna is

A(θ) = -min 12 θ θ 3dB 2
, 30 dB, θ 3dB = 70 • .

We assume that the UTs are divided into G = 2 groups as described in [START_REF] Adhikary | Joint spatial division and multiplexing-the large-scale array regime[END_REF]. The pathloss between UT m in the group g and cell ℓ is selected as in [START_REF] Adhikary | Joint spatial division and multiplexing-the large-scale array regime[END_REF]:

PL(d j,m ) = 1/ 1+ dj,m d0
δ , where δ = 3.7 is the pathless exponent and d 0 = 30 m is the reference distance. We use the channel covariance model in ( 3), but scale it by the pathloss and antenna radiation pattern.

For simplicity, we let the effective pilot SNR, ρ tr , and all the downlink SNRs, ρ dl = P ℓ σ 2 , be 15 dB. We compare the average user throughput,

1 KL L j=1 K m=1 E[log 2 (1 + γ j,m )],
of the proposed TPE precoding with that of conventional RZF precoding. The coefficients of TPE precoding are achieved by solving the relaxed problem in (23) offline. Using Monte-Carlo simulations, Fig. 2 shows the performance for K = 40 users and different number of antennas at each BS: M ∈ {80, 160, 240, 320, 400}. Unlike the single-cell case analyzed in [START_REF] Müller | Linear precoding based on truncated polynomial expansion-Part I: Large-scale single-cell systems[END_REF], where TPE precoding is merely an approximation of RZF precoding, we see that TPE precoding achieves higher user rates for all J j ≥ 3 in the multi-cell case. This can be attributed to the tractable offline optimization of the polynomial coefficients which allows for a better coordination of intra-cell and inter-cell interference, a feature which could not be implemented for the RZF precoding.

CONCLUSION

In this paper, we generalize the low-complexity TPE precoding family from [START_REF] Müller | Linear precoding based on truncated polynomial expansion-Part I: Large-scale single-cell systems[END_REF] to multi-cell scenarios. In particular, we derive deterministic equivalents for the asymptotic SINRs in massive MIMO systems where the number of antennas and users grow large with a fixed ratio. The most interesting feature of these expressions is that they only depend on the channel statistics and not on the instantaneous channel realizations. This enables us to optimize TPE precoding in an offline manner. The performance of the proposed precoding method is illustrated using simulations. We note that contrary to the single-cell case, TPE precoding outperforms the RZF precoding in terms of both complexity and throughput.

  M . The coefficients are optimized later in this paper. The complexity of TPE precoding scales only as M K[START_REF] Müller | Linear precoding based on truncated polynomial expansion-Part I: Large-scale single-cell systems[END_REF].

  ,m are functions of the pth derivatives δ (p) ℓ,k and T (p) ℓ of δ ℓ,k (t) and T ℓ (t) at t = 0, p = 0, . . . , n. Denote by X (0) j,m = 1

Fig. 1 .

 1 Fig. 1. Illustration of the three-sector site deployment with L = 3 cells considered in the simulations.

5 Fig. 2 .

 52 Fig. 2. Comparison between conventional RZF precoding and the proposed TPE precoding with different orders J = J j ∀j. max W1,...,W L ,ξ ξ (23)


N Chauffert
email: nicolas.chauffert@gmail.com

P Ciuciu

P Weiss

Gradient waveform design for variable density sampling in Magnetic Resonance Imaging

Nicolas Chauffert, Pierre Weiss, Jonas Kahn and Philippe Ciuciu.

Abstract-Fast coverage of k-space is a major concern to speed up data acquisition in Magnetic Resonance Imaging (MRI) and limit image distortions due to long echo train durations. The hardware gradient constraints (magnitude, slew rate) must be taken into account to collect a sufficient amount of samples in a minimal amount of time. However, sampling strategies (e.g., Compressed Sensing) and optimal gradient waveform design have been developed separately so far. The major flaw of existing methods is that they do not take the sampling density into account, the latter being central in sampling theory. In particular, methods using optimal control tend to agglutinate samples in high curvature areas. In this paper, we develop an iterative algorithm to project any parameterization of k-space trajectories onto the set of feasible curves that fulfills the gradient constraints. We show that our projection algorithm provides a more efficient alternative than existinf approaches and that it can be a way of reducing acquisition time while maintaining sampling density for piece-wise linear trajectories.

Index Terms-gradient waveform design, k-space trajectories, variable density sampling, gradient hardware constraints, magnetic resonance imaging

I. INTRODUCTION

T HE advent of new hardware and sampling theories pro- vide unprecedented opportunities to reduce acquisition times in MRI. The design of gradient waveforms minimizing the acquisition time while providing enough information to reconstruct distortion-free images is however an important challenge. Ideally, these two concerns (sampling scheme and gradient waveform designs) should be addressed simultaneously but it is unclear how to design a feasible waveform corresponding to a k-space sampling with good space coverage properties. The gradient waveform design is thus generally performed sequentially: a first step aims at finding the trajectory support or at least control points, and a second step consists of designing the fastest gradient waveforms to traverse the trajectory.

Recent progresses in compressed sensing (CS) [START_REF] Puy | On variable density compressive sampling[END_REF][START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF][START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF][START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF] give hints concerning the best strategies of k-space sampling. In particular, the k-space has to be sampled with a variable density, decaying from the lowest to the highest frequencies. To date, CS sampling schemes are limited to simple trajectories, such as 2D independent sampling and acquisition in the orthogonal readout direction. The design of physically plausible acquisition schemes with a prescribed sampling density remains an open question. To better understand how the proposed approach goes beyond the state-of-the-art, we start with a short review of existing methods for designing gradient waveforms for k-space sampling. We show that existing methods do not encompass the case of curves with high curvature in terms of sampling density and scanning time, justifying a new strategy for designing gradient waveforms.

A. Related works Some heuristic methods have been proposed to design simultaneously the sampling density and the gradient waveforms [START_REF] Mir | Fast three-dimensional k-space trajectory design using missile guidance ideas[END_REF][START_REF] Dale | Optimal design of kspace trajectories using a multi-objective genetic algorithm[END_REF] but these are not computationally efficient and it is not clear how to sample the k-space according to a prescribed distribution using these methods. Let us mention other algorithms that have been proposed to traverse specific curve shapes, such as spiral [START_REF] Gurney | Design and analysis of a practical 3D cones trajectory[END_REF][START_REF] Pipe | Spiral trajectory design: A flexible numerical algorithm and base analytical equations[END_REF] or concentric rings [START_REF] Wu | MRI using a concentric rings trajectory[END_REF].

It is in general easier to consider a sampling trajectory in a first time and design the gradient waveform afterwards. To the best of our knowledge, the latter problem is currently solved using convex optimization [START_REF] Simonetti | An optimal design method for magnetic resonance imaging gradient waveforms[END_REF][START_REF] Hargreaves | Time-optimal multidimensional gradient waveform design for rapid imaging[END_REF], optimal control [START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF], or optimal interpolation of k-space control points [START_REF] Davids | Fast and robust design of time-optimal k-space trajectories in MRI[END_REF]. Given an input curve, both methods find a parameterization that minimizes the time required to traverse the curve subject to kinematic constraints. This principle suffers from two limitations in certain situations. First, reparameterizing the curve changes the density of samples along the curve. This density is now known to be a key factor in compressed sensing [START_REF] Puy | On variable density compressive sampling[END_REF][START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF][START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF][START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF], since it directly impacts the number of required measurements to ensure accurate reconstruction. Second, the challenge of rapid acquisitions is to reduce the scanning time (echo train duration) and limit geometric distortions induced by inhomogeneities of the static magnetic field (B 0 ) by covering the kspace as fast as possible. The perfect fit to an arbitrary curve may be time consuming, especially in the high curvature parts of the trajectory. In particular, the time to traverse piecewise linear sampling trajectories [START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF][START_REF] Chauffert | Travelling salesmanbased compressive sampling[END_REF][START_REF] Chauffert | From variable density sampling to continuous sampling using Markov chains[END_REF][START_REF] Wang | Smoothed random-like trajectory for compressed sensing MRI[END_REF][START_REF] Willett | Errata: Sampling trajectories for sparse image recovery[END_REF] may become too long. Indeed, the magnetic field gradients have to be set to zero at each singular point of such trajectory. For this reason, new methods have to be pushed forward to fill this gap.

B. Contributions

In this paper, we propose an alternative method based on a convex optimization formulation. Given any parameterized curve, our algorithm returns the closest curve that fulfills the gradient constraints. The main advantages of the proposed approach are the following: i) the time to traverse the kspace is usually shorter, ii) the distance between the input and output curves is the quantity to be minimized ensuring a low deviation to the original sampling distribution, iii) the maximal acquisition time may be fixed, enabling to find the closest curve in a given time and iv) it is flexible enough to handle additional hardware constraints (e.g., trajectory starting from the k-space center) in the same framework.

C. Paper organization

In Section II, we review the formulation of MRI acquisition, by recalling the gradient constraints and introducing the projection problem. Then, in Section III, it is shown that curves generated by the proposed strategy (initial parameterization + projection onto the set of physical constraints) may be used to design MRI sampling schemes with locally variable densities. In Section IV, we provide an optimization algorithm to solve the projection problem, and estimate its rate of convergence. Finally, we illustrate the behaviour of our algorithm on two particular cases: rosette and Travelling Salesman Problembased trajectories (Section V).

II. DESIGN OF k-SPACE TRAJECTORIES USING PHYSICAL

GRADIENT WAVEFORMS.

In this section, we recall classical modelling of the acquisition constraints in MRI [START_REF] Hargreaves | Time-optimal multidimensional gradient waveform design for rapid imaging[END_REF][START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF]. We justify the lack of accuracy of current reparameterization methods in the context of variable density sampling, and motivate the introduction of a new projection algorithm.

A. Sampling in MRI

In MRI, images are sampled in the k-space domain along parameterized curves s : [0, T ] → R d where d ∈ {2, 3} denotes the image dimensions. The i-th coordinate of s is denoted s i . Let u : R d → C denote a d dimensional image and û be its Fourier transform. Given an image u, a curve s : [0, T ] → R d and a sampling step ∆t, the image u shall be reconstructed using the set 1 :

E = û(s(j∆t)), 0 j T ∆t . (1)

B. Gradient constraints

The gradient waveform associated with a curve s is defined by g(t) = γ -1 ṡ(t), where γ denotes the gyro-magnetic ratio [START_REF] Hargreaves | Time-optimal multidimensional gradient waveform design for rapid imaging[END_REF]. The gradient waveform is obtained by energizing gradient coils (arrangements of wire) with electric currents.

1) kinematic constraints: due to obvious physical constraints, these electric currents have a bounded amplitude and cannot vary too rapidly (slew rate). Mathematically, these constraints read:

g G max and ġ S max
where • denotes either the ∞ -norm defined by f ∞ := max 1≤i≤d sup t∈[0,T ] |f i (t)|, or the ∞,2 -norm defined by 1 For ease of presentation, we assume that the values of u in the k-space correspond to its Fourier transform and we neglect distorsions occuring in MRI such as noise.

f ∞,2 := sup t∈[0,T ] d i=1 |f i (t)| 2 1 2 . These constraints might be Rotation Invariant (RIV) if • = • ∞,2 or Rotation Variant (RV) if • = • ∞ ,
depending on whether each gradient coil is energized independently from others or not. The set of kinematic constraints is denoted S:

S := s ∈ C 2 ([0, T ]) d , ṡ α, s β . (2)
2) Additional affine constraints: Specific MRI acquisitions may require additional constraints, such as:

• Imposing that the trajectory starts from the k-space center (i.e., s(0) = 0) to save time and avoid blips. The endpoint can also be specified by s(T ) = s T , if s T can be reached during travel time T . • In the context of multi-shot MRI acquisition, several radio-frequency pulses are necessary to cover the whole k-space. Hence, it makes sense to enforce the trajectory to start from the k-space center at every T R (repetition time):

s(k • T R) = 0, 0 k T T R .
• In addition to starting from the k-space center, one could impose starting at speed 0, i.e. ṡ(0) = 0. • To avoid artifacts due to flow motion in the object of interest, gradient moment nulling (GMN) techniques have been introduced in [START_REF] Majewski | MRI pulse sequence design with first-order gradient moment nulling in arbitrary directions by solving a polynomial program[END_REF]. In terms of constraints, nulling the i th moment reads t t i g(t)dt = 0. Each of these constraints can be modeled by an affine relationship. Hereafter, the set of affine constraints is denoted by A:

A := s : [0, T ] → R d , A(s) = v ,
where v is a vector of parameters in R p (p is the number of additional constraints) and A is a linear mapping from the curves space to R p . A sampling trajectory s : [0, T ] → R d will be said to be admissible if it belongs to the set S ∩ A. In what follows, we assume that this set is non-empty, i.e. S ∩A = ∅. Moreover, we assume, without loss of generality, that the linear constraints are independent (otherwise some could be removed).

C. Finding an optimal reparameterization

The traditional approach to design an admissible curve s ∈ S given an arbitrary curve c : [0, T ] → R d consists of finding a reparameterization p such that s = c • p satisfies the physical constraints while minimizing the acquisition time. This problem can be cast as follows:

T Rep = min T such that ∃ p : [0, T ] → [0, T ], c • p ∈ S. (3)
It can be solved efficiently using optimal control [START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF] or convex optimization [START_REF] Hargreaves | Time-optimal multidimensional gradient waveform design for rapid imaging[END_REF]. The resulting solution s = c • p has the same support as c. This method however suffers from an important drawback when used in the CS framework: it does not provide any control on the density of samples along the curve. For example, for a given curve support shown in Fig. 1(a), we illustrate the new parameterization (keeping the same support) and the corresponding magnetic field gradients (see Fig. 1(b) for a discretization of the curve and (c) for the gradient). We notice that the new parameterized curve has to stop at every angular point of the trajectory, yielding more time spent by the curve in the neighbourhood of these points (and more points in the discretization of the curve in Fig. 1(b)). This phenomenon is likely to modify the sampling distribution, as illustrated in Section III.

The next part is dedicated to introducing an alternative method relaxing the constraint of keeping the same support as c.

D. Projection onto the set of constraints

The idea we propose in this paper is to find the projection of the given input curve c onto the set of admissible curves S:

s * := argmin s∈S 1 2 d 2 (s, c) = argmin s∈S 1 2 s -c 2 2 ( 4
)
where

d 2 (s, c) = s -c 2 2 := T t=0 s(t) -c(t) 2 2 dt
. This method presents important differences compared to the optimal control approach introduced here above: i) the solution s * and c have different support (see Fig. 1(d)) unless c is admissible; ii) the sets composed of the discretization of c and s * at a given sampling rate are close to each other (Fig. 1(e)); iii) the acquisition time T is fixed and equal to that of the input curve c. In particular, time to traverse a curve is in general shorter than using exact parameterization (see Fig. 1(f) where

T < T Rep ). (a) (b) (c) g(t) t (d) (e) (f) g(t)
t In the next section, we will explain why the empirical distribution of the curves obtained by this projection method is close to the one of the input curve, and we illustrate how the parameterization will distort the sampling distribution.

III. CONTROL OF THE SAMPLING DENSITY

Recent works have emphasized the importance of the sampling density [START_REF] Puy | On variable density compressive sampling[END_REF][START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF][START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF][START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF] in an attempt to reduce the amount of acquired data while preserving image quality at the reconstruction step. The choice of an accurate distribution p is crucial since it directly impacts the number of measurements required [START_REF] Rauhut | Compressive Sensing and Structured Random Matrices[END_REF]. In this paper, we will denote by π a distribution defined over the k-space K. The profile of this distribution can be obtained by theoretical arguments [START_REF] Puy | On variable density compressive sampling[END_REF][START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF][START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF][START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF] leading to distributions as the one depicted in Fig. 2(a). Some heuristical distributions are known to perform well in CS-MRI experiments (Fig. 2(b)). A comparison between these two approaches can be found in [START_REF] Chauffert | Variable density compressed sensing in MRI. Theoretical vs. heuristic sampling strategies[END_REF]. However, designing a trajectory that performs sampling according to a fixed distribution while satisfying gradient constraints, is really challenging and has never been adressed so far. Recent attemps did not manage gradient constraints to design continuous variable density sampling trajectories since the sampling curve cannot be traversed at constant speed, especially at angular points as shown in Fig. 1(b). This impacts the sampling density by concentrating probability mass in the portions of the curve associated with large curvature. In this section, we show that if we start from an inadmissible trajectory that covers the k-space with an empirical distribution close to the target one, then we can figure out how to design an admissible sampling trajectory with almost the same empirical distribution. Strategies to design such input curves are detailed in Appendix A. In what follows, we first derive theoretical guarantees, and then, we perform experiments to quantify the distribution distortion.

A. Theoretical guarantees

Let us show that the sampling distribution carried over a curve that our algorithm delivers can be close to a fixed target distribution. To this end, we start with the definition of the empirical distribution of a curve, and then we introduce a distance defined over distributions.

Definition 1 (Empirical measure of a curve). Let γ denote the Lebesgue measure and γ T = γ T denote the Lebesgue measure normalized on the interval [0, T ]. The empirical measure of a curve s : [0, T ] → K ⊆ R d is defined for any measurable set ω of K as:

P s (ω) = γ T (s -1 (ω)).
This means that the mass of a set ω is proportional to the time spent by the curve in ω. The natural distance arising in our work is the between-distribution Wasserstein distance W 2 defined hereafter: Definition 2 (Wasserstein distance W 2 ). Let M be a domain of R d and P(M ) be the set of measures over M . For µ, ν ∈ P(M ), W 2 is defined as:

W 2 (µ, ν) = inf γ∈Π(µ,ν)
x -y 2 2 dγ(x, y)

(5)

where Π ⊂ P(M × M ) denote the set of measures over M × M with marginals µ and ν on the first and second factors, respectively.

W 2 is a distance over P(M ) (see e.g., [START_REF] Villani | Optimal transport: old and new[END_REF]). Intuitively, if µ and ν are seen as mountains, the distance is the minimum cost of moving the mountains of µ into the mountains of ν, where the cost is the 2 -distance of transportation multiplied by the mass moved. Hence, the coupling encodes the deformation map to turn one distribution into the other.

A global strategy to design feasible trajectories with an empirical distribution close to a target density π relies on:

1) finding an input curve c, such that P c is close to π; 2) estimating s * the projection of c onto the set of constraints, by solving Eq. ( 4). The objective is to show that W 2 (P s * , π) is small. Since W 2 is a norm, the triangle inequality holds:

W 2 (P s * , π) W 2 (P s * , P c ) Distortion + W 2 (P c , π) Initial guess . ( 6
)
The initial guess term W 2 (P c , π) can be as small as possible if c is a Variable Density Sampler [START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF]. The proof is postponed to Appendix A for the ease of reading. Here, we are interested in the distortion term W 2 (P s * , P c ). The following proposition shows that the W 2 distance between the empirical distributions of the input curve and of the output curve (c and s * ) is controlled by the quantity d(s * , c) to be minimized when solving Eq. ( 4). Proof. In terms of distributions, the quantity d(s, c) reads:

d 2 (s, c) = M ×M x -y 2 2 dγ s,c (x, y) (7)
where γ s,c is the coupling between the empirical measures P s and P c defined for all couples of measure sets

(ω 1 , ω 2 ) ∈ M 2 by γ c,s (ω 1 , ω 2 ) = 1 T T t=0 1 ω1 (s(t))1 ω2 (c(t))dt.
The choice of this coupling is equivalent to choosing the transformation map as the association of locations of c(t) and s(t) for each t. We notice that the quantity to be minimized in Eq. ( 7) is an upper bound of W 2 (P s , P c ) 2 , with the specific coupling γ s,c .

To sum up, solving the projection problem amounts to minimizing an upper bound of the W 2 distance between the target and the empirical distributions of the solution if we neglect the influence of the initial guess term.

B. Simulations

Next, we performed simulations to show that the sampling density is better preserved using our approach than relying on the optimal control approach. For doing so, we use travelling salesman-based (TSP) sampling trajectories [START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF][START_REF] Chauffert | From variable density sampling to continuous sampling using Markov chains[END_REF], which are an original way to design random trajectories which empirical distribution is any target distribution such as in Fig. 3(a). 10000 such independent TSP were drawn and parameterized with arc-length: note that these parameterizations are not admissible in general. Then, we sampled each trajectory at fixed rate ∆t (as in Fig. 3(b)), to form an histogram (the empirical distribution shown in Fig. 3(c)), which was compared to π in Fig. 3(d). The error is actually not closed to zero, since the convergence result is asymptotic in the size of the curve that is bounded in this experiment. In Fig. 4 (top row), we show that the classical reparameterization technique [START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF] lead to a major distortion of the sampling density, because of the behavior on the angular points illustrated in Fig. 1(b). Then, we considered three constant speed parameterizations and projected them onto the same set of constraints (G max = 40 mT.m -1 and S max = 150 mT.m -1 .ms -1 ). Among these 3 initial guesses, we first used an initial parameterization with low velocity (10 % of the maximal speed γG max with γ = 42.576 MHz.T -1 ), which projection fits the sampling density quite well. Then, we increased the velocity to reach first 50 % and then 100 % of the maximal speed. The distortion of the sampling density of the projected curve increased, but remained negligible in contrast to the exact reparameterization. Hence, this example illustrates that if we have a continuous trajectory with an empirical sampling distribution close to the target one, our projection algorithm enables to design feasible waveforms while sampling the k-space along a discretized trajectory with an empirical density close to the target one as well.

IV. FINDING FEASIBLE WAVEFORMS USING CONVEX

OPTIMIZATION

Since the set of constraints S ∩ A is convex, closed and non-empty, Problem (4) always admits a unique solution. Even though S has a rather simple structure 2 , it is unlikely that an explicit solution to Problem (4) can be found. In what follows, we thus propose a numerical algorithm to find the projection.

Problem discretization: In order to define a numerical algorithm, we first discretize the problem. A discrete curve s is defined as a vector in R n•d where n is the number of discretization points. Let s(i) ∈ R d denote the curve location at time (i -1)δt with δt = T n-1 . The discrete derivative ṡ ∈ R n•d is defined using first order differences:

ṡ(i) = 0 if i = 1, (s(i) -s(i -1))/δt if i ∈ {2, . . . , n}.
In the discrete setting, the first-order differential operator can be represented by a matrix Ṁ ∈ R nd×nd , i.e. ṡ = Ṁs. We define the discrete second-order differential by M = -Ṁ * Ṁ ∈ R nd×nd .

2 it is just a polytope when the ∞ -norm is used.

An efficient projection algorithm: The discrete problem we consider is the same as problem ( 4) except that all objects are discretized. The set S is now S := {s ∈ R n•d , Ṁs α, Ms β} with all norms discretized. Similarly, the discretized version of A is denoted A. The main idea of our algorithm is to take advantage of the structure of the dual problem to design an efficient projection algorithm. The following proposition specifies the dual problem and the primal-dual relationships.

Proposition 2. Let q * := sup q ≤1 q, q denote the dual norm of • . The following equality holds

min s∈S∩A 1 2 s -c 2 2 (8) = sup q1,q2∈R n•d F (q 1 , q 2 ) -α q 1 * -β q 2 * , (9)
where

F (q 1 , q 2 ) = min s∈A Ṁs, q 1 + Ms, q 2 + 1 2 s -c 2 2 . ( 10
)
Moreover, let (q * 1 , q * 2 ) denote any minimizer of the dual problem (9), s * denote the unique solution of the primal problem (8) and s * (q * 1 , q * 2 ) denote the solution of the minimization problem [START_REF] Simonetti | An optimal design method for magnetic resonance imaging gradient waveforms[END_REF]. Then s * = s * (q * 1 , q * 2 ). The following proposition gives an explicit expression of s * (q * 1 , q * 2 ). Proposition 3. The minimizer s * (q * 1 , q * 2 ) = arg min s∈A Ṁs, q 1 + Ms, q 2 + 1 2 sc 2 2 is given by

s * (q 1 , q 2 ) = z + A + (v -Az), ( 11
)
where

A + = A * (AA * ) -1 denotes the pseudo-inverse of A and z = c -Ṁ * q 1 -M * q 2 .
Let us now analyze the smoothness properties of F .

Proposition 4. Function F (q 1 , q 2 ) is concave differentiable with gradient given by

∇F (q 1 , q 2 ) = - Ṁs * (q 1 , q 2 ) Ms * (q 1 , q 2 ) . ( 12
)
Moreover, the gradient mapping ∇F is Lipschitz continuous with constant L = ||| Ṁ * Ṁ+ M * M|||, where |||M||| denotes the spectral norm of M.

The proofs are given in Appendix B. The dual problem (9) has a favorable structure for optimization: it is the sum of a differentiable convex function F (q 1 , q 2 ) = -F (q 1 , q 2 ) and of a simple convex function G(q 1 , q 2 ) = α q 1 * +β q 2 * . The sum F+G can thus be minimized efficiently using accelerated proximal gradient descents [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF] (see Algorithm 1 below).

Moreover, by combining the convergence rate results of [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF][START_REF] Beck | Gradient-based algorithms with applications to signal recovery[END_REF] and some convex analysis (see Appendix D), we obtain the following convergence rate: Algorithm 1: Projection algorithm in the dual space

Input: c ∈ R n•d , α, β > 0, n it . Output: s ∈ R n•d an approximation of the solution s * . Initialize q (0) = (q (0) 1 , q (0) 2 ) with q (0) i = 0 for i = 1, 2. Set y (0) = q (0) . Set t = 1/L. for k = 1 . . . n it do q (k) = prox tG (y (k-1) -t∇ F (y (k-1) )) y (k) = q (k) + k-1 k+2 (q (k) -q (k-1) ) return s = s * q (nit) 1
, q (nit) 2 ; Theorem 1. Algorithm 1 ensures that the distance to the minimizer decreases as O 1 k 2 :

s (k) -s * 2 2 ≤ 2L q (0) -q * 2 2 k 2 . (13)
V. ILLUSTRATIONS To compare our results with [START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF], we used the same gradient constraints. In particular, the maximal gradient norm G max set to 40 mT.m -1 , and the slew-rate S max to 150 mT.m -1 .ms -1 . k-space size is 6 cm -1 and sampling rate is ∆t = 4 µs. For the ease of trajectory representation, we limit ourselves to 2D sampling curves, but the algorithm encompasses the 3D setting. The Matlab codes containing the projection algorithm as well as the scripts to reproduce the results depicted hereafter are available at http://chauffertn.free.fr/codes.html. First, we illustrate the output of Algorithm 1 for the case of rosette trajectories introduced in [START_REF] Noll | Multishot rosette trajectories for spectrally selective MR imaging[END_REF], showing a similar behaviour of our algorithm compared to optimal control. Then, we show an application to gradient waveform design to traverse TSP-based sampling trajectories. In this case, only our algorithm yields a fast kspace coverage.

A. Smooth trajectories

The first illustration in Fig. 5 shows that when the sampling curve is smooth enough, our algorithm has a similar behaviour than optimal control. Here, we considered the case of Rosette trajectory defined in [START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF][START_REF] Noll | Multishot rosette trajectories for spectrally selective MR imaging[END_REF] as k(t) = k max sin(ω 1 t)[cos(ω 2 t); sin(ω 2 t)] (with ω 1 = 1.419 and ω 2 = 0.8233 rad.s -1 ). In Fig. 5, we considered c(t) the parameterization of k(t) at constant speed (90% of γG max ). The black curve represents the projection of c onto the set of constraints for Rotation Variant (RV) constraints (Fig. 5(a)) and for Rotation Invariant (RIV) constraints (Fig. 5(b)). Then, we show the corresponding gradient waveforms to traverse the trajectory (see Fig. 5(c,d) for optimal control solution and Fig. 5(e,f) for the output of our projection algorithm).

We observed that the time spent to traverse the curve is slighty shorter compared to the exact reparameterization, while the support of the curve is reduced. Of course, the initial parameterization of c has a strong impact, and there is a trade-off between the time spent to traverse the curve and the distance to the original support. We also noticed that RV constraints are less restrictive, enabling a faster traversal time in the optimal control framework (Fig. 5(c)) or a lower support distortion in the proposed setting (Fig. 5(a)), as also noticed in [START_REF] Vaziri | The fastest gradient waveforms for arbitrary and optimized k-space trajectories[END_REF].

In this example, there is no real difference between the optimal control and the projection approaches. The following example, based on piece-wise linear sampling curves, shows in contrast that optimal control yields really longer traversal time. Fig. 5. Gradient waveforms to design rosette trajectories [START_REF] Noll | Multishot rosette trajectories for spectrally selective MR imaging[END_REF]. Left (resp. right), with Rotation-Variant (RV) (resp. Rotation-Invariant (RIV)) constraints ( = ∞ (resp. = ∞,2 )). Corresponding gradient waveforms (gx(t), gy(t), g(t) ). Top row: k-space trajectories between -6 and 6 cm -1 . Middle row: gradient waveforms corresponding to optimal reparameterization [START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF]. Bottom row: gradient waveforms corresponding to the proposed projection algorithm.

B. TSP sampling

Next, we performed similar experiments with a TSP trajectory, as introduced in [START_REF] Chauffert | From variable density sampling to continuous sampling using Markov chains[END_REF]. We illustrate the results for RIV constraints, given an initial parameterization at 50 % of the maximal speed. In Fig. 6(a), we notice that the output trajectory is a smoothed version of the initial piece-wise linear trajectory. Algthough the output curve has a slightly different support, the traversal time is reduced from 58 ms (optimal reparameterization that can be computed exactly) to 16 ms (Fig. 6(b-c)). The main reason is that TSP trajectories embody singular points, requiring the gradients to be set to zero for each of these points. Therefore, a sampling trajectory with singular points is time consuming. The main advantage of our algorithm is that the trajectory can be smoothed around these points, which saves a lot of acquisition time.

This second example shows that with existing methods, it is hopeless to implement TSP-based sequences in many MRI modalities, since the time to collect data can be larger than any realistic repetition time. In contrast, our method enables to traverse such curves in a reasonable time which can be tuned according to the image weighted (T 1 , T 2 , proton density). (T = 16 ms) Fig. 6. Gradient waveforms to design TSP trajectories [START_REF] Noll | Multishot rosette trajectories for spectrally selective MR imaging[END_REF]. (a): k-space trajectories between -6 and 6 cm -1 . (b): gradient waveforms (g = (gx, gy)) corresponding to optimal reparameterization [START_REF] Lustig | A fast method for designing timeoptimal gradient waveforms for arbitrary k-space trajectories[END_REF]. (c): gradient waveforms corresponding to the proposed projection algorithm.

VI. CONCLUSION

We developed an algorithm to project any parameterized curve onto the set of curves which can be implemented on actual MRI scanners. Our method is an alternative to the existing gradient waveform design based on optimal control. The main advantage of our approach is that one can reduce the time to traverse a curve if the latter has a lot of high curvature areas, like a solution of the TSP.

A. Definition of a VDS First, we need to introduce the definition of weak convergence for measure: Definition 3. A sequence of measures µ n ∈ P(K), the set of distributions defined over K, is said to weakly converge to µ if for any bounded continuous function φ

K φ(x)dµ n (x) → K φ(x)dµ(x).
We use the notation µ n µ.

According to [START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF], a (generalized) π-VDS is a set of times T n , such that T n → ∞ when n → ∞, and a sequence of curves c Tn : [0, T n ] → R d such that P c Tn π when n tends to infinity. A consequence of the definition is that the relative time spent by the curve in a part of the k-space is proportional to its density. Before showing that this implies that W 2 (P c Tn , π) tends to 0, we give two examples of VDS.

B. VDS examples

We give two examples to design continuous sampling trajectories that match a given distribution. The two examples we propose provide a sequence of curves, hence a sequence of empirical measures, that weakly converge to the target density.

1) Spiral sampling: The spiral-based variable density sampling is now classical in MRI [START_REF] Spielman | Magnetic resonance fluoroscopy using spirals with variable sampling densities[END_REF]. For example, let n ∈ N be the number of revolutions and r : [0, 1] → R + a strictly increasing smooth function. Denote by r -1 its inverse function.

Define the spiral for t ∈ [0, n] by c n (t) = r t n cos(2πt) sin(2πt) and the target distribution π by:

π(x, y) =      . r -1 √ x 2 +y 2 2π r(1) ρ=r(0) . r -1 (ρ)ρdρ if r(0)
x 2 +y 2 r(1) 0 otherwise then P cn π when the number of revolutions n tends to infinity.

2) Travelling Salesman-based sampling: The idea of using the shortest path amongst a set of points (the "cities") to design continuous trajectories with variable densities has been justified in [START_REF] Chauffert | Variable density sampling with continuous trajectories. Application to MRI[END_REF][START_REF] Chauffert | Travelling salesmanbased compressive sampling[END_REF]. Let us draw n k-space locations uniformly according to a density q define over the dD k-space (d = 2 or 3), and join them by the shortest path (the Travelling Salesman solution). Then, denote by c n a constant-speed parameterization of this curve. Define the density:

π = q (d-1)/d q (d-1)/d (x)d(x)
Then P cn π when the number of cities n tends to infinity. These two sampling strategies are efficient to cover the kspace according to target distributions, as depicted in Fig. 7. For spiral sampling, the target distribution may be any 2D radial distribution, whereas the Travelling salesman-based sampling enable us to consider any 2D or 3D density.

C. Control of W 2 distance

Let us now assume without loss of generality that K =

[-k max , k max ] d .
Let us recall a central result about W 2 (see e.g., [START_REF] Villani | Optimal transport: old and new[END_REF]):

Proposition 5. Let M ⊂ R d , µ ∈ P(M ) and µ n be a sequence of P(M ). Then, if M is compact µ n µ ⇔ W 2 (µ n , µ) → 0
An immediate consequence of this proposition and of the compactness of K is the following proposition: Proposition 6. Let (c Tn ) n 1 be a π-VDS, and ε > 0. Then, there exists n 1 such that c Tn : [0, T n ] → K fulfills:

W 2 (P c Tn , π) ε.
To sum up, Proposition 6 ensures that we can find an input curve which empirical distribution is as close to the target distribution π as we want. Ṁs, q 1 -α q 1 * + Ms, q 2 -β q 2 = sup q1,q2∈R -α q 1 * -β q 2 *

The relationship between the primal and dual solutions reads s * = arg min s∈A 1 2 sc 2 2 + s, Ṁ * q * 1 + s, M * q * 2 . The sup and the min can be interverted at the third line, due to standard theorems in convex analysis (see e.g. [START_REF] Rockafellar | Convex analysis[END_REF]Theorem 31.3]).

APPENDIX C PROOF OF PROPOSITIONS 4 AND 3

To show Proposition 3, first remark that arg min s∈A Ṁs, q 1 + Ms, q 2 + 1 2 sc 2 2 = arg min

s∈A 1 2 s -(c -Ṁq 1 -M * q 2 ) 2 2 .
Therefore, s * (q 1 , q 2 ) is the orthogonal projection of z = c -Ṁq 1 -M * q 2 onto A. Since A is not empty, AA + v = v, and the set A = {s ∈ R n•d , As = v} can be rewritten as

A = A + v + ker(A).
The vector zs * (q 1 , q 2 ) is orthogonal to A, it therefore belongs to ker(A) ⊥ = im(A * ). Thus s * (q 1 , q 2 ) = z + A * λ for some λ such that:

A(z + A * λ) = v.
This leads to λ = (AA * ) -1 (v -Az). We finally get s * (q 1 , q 2 ) = z + A * (AA * ) -1 (v -Az), ending the proof. Proposition 4 results from standard results of convex optimization. See for instance [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]Theorem 1] or the book [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms I: Part 1: Fundamentals[END_REF].

APPENDIX D PROOF OF THEOREM 1.

The analysis proposed here follows closely ideas proposed in [START_REF] Beck | Gradient-based algorithms with applications to signal recovery[END_REF][START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF][START_REF] Boyer | An algorithm for variable density sampling with block-constrained acquisition[END_REF][START_REF] Beck | A fast dual proximal gradient algorithm for convex minimization and applications[END_REF]. We will need two results. The first one is a duality result from [START_REF] Boyer | An algorithm for variable density sampling with block-constrained acquisition[END_REF]. Proposition 8. Let f : R m → R ∪ {∞} and g : R n → R ∪ {∞} denote two closed convex functions. Assume that g is σstrongly convex [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms I: Part 1: Fundamentals[END_REF] and that Ari(dom(f ))∩ri(dom(g)) = ∅. Let A ∈ R m×n denote a matrix.

Let p(x) = f (Ax) + g(x) and d(y) = -g * (A * y) -f * (y). Let x * be the unique minimizer of p and y * be any minimizer of d.

Then g * is differentiable with 1 σ Lipschitz-continuous gradient. Moreover, by letting x(y) = ∇g * (-A * y):

x(y) -x * The second ingredient is the standard convergence rate for accelerated proximal gradient descents given in [START_REF] Beck | Gradient-based algorithms with applications to signal recovery[END_REF]. Proposition 9. Under the same assumptions as Proposition 8, consider Algorithm 2.

Then

y (nit) -y * 2 2 = O |||A||| 2 σ•n 2 it .
To conclude, it suffices to set g(s) = 1 2 s-c 2 2 , f (q 1 , q 2 ) = ı Bα (q 1 ) + ı Bα (q 2 ) and A = Ṁ M . By doing so, the

Fig. 1 .

 1 Fig. 1. Comparison of two methods to design gradient waveforms. Top row: Optimal control-based parameterization [12]. (a): input curve support. (b): discrete representation of the optimal reparameterization of the curve in S. (c): corresponding gradient waveforms (gx, gy). Dashed lines correspond to 0 and +/-Gmax. Bottom row: Illustration of the projection algorithm. (d): the same input curve as in (a) parameterized at maximal speed, and the support of the projected curve onto the set S. (e): discrete representation of the input and projected curves. (f): corresponding gradient waveforms (gx, gy) with the same time scale as in (c): time to traverse the curve is reduced by 39%.

Fig. 2 .

 2 Fig. 2. Examples of 2D sampling distribution. (a): optimal distribution for a Symmlet transform [1, 2]. (b): radial distribution advocated in [21]: p(k) ∝ (1 -|k|/kmax) p with p = 3.

Proposition 1 .

 1 For any two curves s and c : [0, T ] → R d : W 2 (P c , P s ) d(s, c).

Fig. 3 .

 3 Fig. 3. Illustration of TSP VDS. (a) Target distribution. (b) Regular sampling of TSP curve parameterized with its arc-length. (c) Empirical distribution (histogram of samples). (d) Histogram difference (c)-(a).

Fig. 4 .

 4 Fig. 4. Illustration of TSP trajectories traversed with optimal control (top row) and with our projection algorithm (rows 2-4). Columns represent the k-space trajectory, the empirical distribution and the difference with the target distribution π (Fig. 3(a)).

Fig. 7 .

 7 Fig. 7. Examples of Variable Density Sampling. (a): TSP trajectory with empirical distribution depicted in Fig. 2(a).(b): spiral trajectory with empirical radial distribution such as Fig. 2(b).

-c 2 2 + 2 s -c 2 2

 222 Definition 4 (indicator function). Let B ⊆ R n . The indicator of B is denoted ı B and defined by:ı B (x) = 0 if x ∈ B +∞ otherwiseLet us now recall a classical result of convex optimization[28, P. 195]:Proposition 7. Let B α = {x ∈ R n , xα}. Then the following identity holds:ı Bα (x) = supy∈R n x, y -α y * . Now, we can prove Proposition 2. ı Bα ( Ṁs) + ı B β ( Ms) = min s∈A 1 + sup q1,q2∈R n•d

  y) -d(y * )).

ACKNOWLEDGMENTS

This work benefited from the support of the FMJH Program Gaspard Monge in optimization and operation research, and from the support to this program from EDF. APPENDIX A DENSITY DEVIATION, CONTROL OF W 2 -DISTANCE.

In Section III-A, we aim at controlling the Wasserstein distance W 2 (P s * , π), where π is a target fixed sampling distribution, and P s * is the empirical distribution of the projected curve. We used the triangle inequlity (6) to bound this quantity by W 2 (P s * , P c ) + W 2 (P c , π). Here, we show that the quantity W 2 (P c , π) can be as small as possible if c is Variable Density Sampler (VDS) [2]. First, we define the concept of VDS, and then we provide two examples. Next, we show that if c is a VDS, W 2 (P c , π) tends to 0 as the length of c tends to infinity.

Algorithm 2: Accelerated proximal gradient descent Input:

projection problem rewrites min s∈R n ḋ p(s) = f (As) + g(s). Its dual problem (9) can be rewritten more compactly as min q=(q1,q2)∈R n ḋ ×R n ḋ d(q) = g * (-A * q) + f * (q). Note that function g is 1-strongly convex. Therefore, Algorithm 2 can be used to minimize d, ensuring a convergence rate in O 1