Pendant ces trois années, j'ai eu l'occasion unique de toucher au monde de la recherche en mathématiques appliquées aux problèmes de la physique et d'améliorer considérablement mes compétences professionnelles. Cela n'aurait pas été possible sans mes directeurs de recherche, les professeurs Paul Armand en mathématiques appliquées et Vincent Kermene, directeur de recherche au département de photonique, qui sont tous deux des professionnels incontestables dans leur domaine.

Le professeur Paul Armand est, sans exagération, un expert de haut niveau. C'est la seule personne que j'ai rencontrée dans ma vie qui peut présenter ses idées d'une manière aussi méthodique, claire et simple et, ce faisant, transmettre la richesse des connaissances qu'il possède à des étudiants comme moi. Dire que je lui suis reconnaissant, c'est ne rien dire. Les qualités personnelles de cet homme sont difficiles à surestimer.

Vincent Kermene, directeur de recherche, est sans aucun doute un professionnel hautement qualifié, avec qui il a été très intéressant et instructif de travailler. Au cours de cette période, mon amour pour la physique et pour l'apprentissage du fonctionnement de notre monde s'est non seulement développé, mais j'ai également acquis de nombreuses connaissances en traitant avec lui. Les qualités personnelles de Vincent, telles que l'ouverture et la convivialité, ont certainement facilité mon adaptation à l'équipe, me permettant d'interagir efficacement avec ses membres et d'obtenir des résultats importants.

Je suis également reconnaissant au reste de l'équipe, à savoir Alain Barthélémy, directeur de recherche émérite, Alexandre Boju, Professeur Agnès Desfarges-Berthelemot, et Geoffrey Maulion, qui appartiennent au département de photonique. Ce fut un plaisir de travailler avec eux pendant cette période. Sans aucun doute, chacun d'entre eux a contribué aux objectifs définis et c'est grâce au travail d'équipe que nous avons pu atteindre nos objectifs.

Cette thèse a été

Chapter 1 Introduction

This work studies algorithms for the dynamic phase control of an array of laser beams [START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF]. The main interest from the physical point of view is to obtain a synthetic beam of high brightness (high power and high beam quality). There exist many physical problems where a beam of high power is required. The first application comes from the military goal to efficiently hit potentially dangerous objects such as self-guided missiles or drones over big distances and is developing in the context of the Tactical Advanced Laser Optical System project (TALOS). The main requirement is that the propagated light must have a high destructive effect, which can be transmitted over a long distance without tangible energy losses. Another application uses a pressure of photons to move light-driven nano-spacecraft in the vacuum of space, which is equipped with a light sail (Breakthrough Starshot), where the same approach can be applied to clean the earth's orbital space from debris (CleanSpace). One more interesting application raises from the need to accelerate particles using a high power laser (ICAN), where the main user is Conseil Européen pour la Recherche Nucléaire (CERN). To achieve the defined task the coherent beam combining [START_REF] Brignon | Coherent Laser Beam Combining[END_REF] approach is applied. The physical fact that is extensively used here is the interference of light waves, which is controlled by their phases. Thus, clearly, if several parallel waves are in phase, in the far field they overlap and exhibit an intense central lobe. An example of such a process is depicted on Figure 1.2, where the near field on the left (an intensity pattern of the four output beams of the laser system) and the far field on the right (an intensity pattern after propagation) are presented. It can be seen that depending on beams phases a coherent in phase beams (high peak at the center) or incoherent, without stable phase relationships (smooth intensity distribution) intensity pattern can be obtained. This problem can be formulated in terms of a mathematical model. Each j-th beam is a light wave with its own amplitude and phase. This information can be encoded into a complex Figure 1.2: Coherent and Incoherent combining example: left, intensity pattern of 4 beams (near field), right, intensity pattern of the their far field when beams are coherent (in phase case) or incoherent (without stable phase relationships) number x j ∈ C, which has an amplitude-phase representation x j = α j • e iϕ j , where α j ≥ 0 equals to a square root of beam intensity, and ϕ j ∈ [-π, π[ represents its phase.

Then, the problem of maximization of brightness for an array of n laser beams can be formulated in terms of the following mathematical model Clearly, it reaches maximum when ϕ 1 = • • • = ϕ n , which coincides with physical interpretation. However, despite the simplicity of the formulation, there are several difficulties, which increase the complexity of the problem.

The first difficulty is that direct phase measurements ϕ 1 , . . . , ϕ n are not possible, which means that these values are not observable directly. The only available possibility is phase modulations by given values δϕ 1 , . . . , δϕ n . One option to receive information about the phases of an array of laser beams is to measure the intensities of the transmitted laser rays through an optical diffuser. The transformation returns the phase-less information which can be used for a vector of phase corrections ϕ ∈ [-π, π[ n computation. This idea is a basis for an opto-numerical algorithm developed in [START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF]. In practice, the algorithm is applied continuously and never stops. The reason comes from the fact that the phases of the beam array change continuously at high frequency up to some kilohertz, which is the case due to the intrinsic noise of the independent amplifiers and because of the sensitivity of the laser system to environmental perturbations such as pressure, temperature or acoustics. Thus, the optical part of the algorithm is done by means of noisy intensity measurements from a photonic system plus phase modulations, and the numerical part is an algorithm, which computes ϕ ∈ [-π, π[ n .

The second difficulty is that we have a need to correct unknown phases ϕ 1 , . . . , ϕ n not just to a constant value to obtain a co-phased far field pattern, but to an arbitrary vector of targets ϕ t 1 , . . . , ϕ t n . This vector describes the atmospheric perturbations which frequently changes and must be taken into account to maintain a high power density far away from the laser system. Thus, it is extremely important for the numerical algorithm to be adaptive and to not require extra computational operations to change the target.

The third difficulty is a presence of a high level of noise in the experimental process. There are at least two possible sources which are the noise during measurements and phase corrections. Thus, it is extremely important that the algorithm, which computes a phase correction vector δϕ, is highly robust.

The crucial point, which is the main topic of this work, is to select from existing numerical algorithms or to develop a new one, that can perform a phase reconstruction satisfying the following criteria:

• Relevancy: the ability of the algorithm to compute a sequence of corrections capable to hit the target signal with high precision.

• Rapidity: the ability to do it in a short time (less than one millisecond).

1.1. NOTATION

• Robustness: the algorithm should be as insensitive as possible to a high level of experimental noises.

• Scalability: the algorithm must remain efficient with an increasing number of input sources.

• Adaptability: dynamic updates of the target signal must be taken into account. The correction process must be able to quickly adapt to a change in the target signal.

Notation

All vector inequalities are understood componentwise. Let x be a vector with complex components. The j-th component of x is denoted by x j . The vectors |x| and arg(x) are the vectors those components are the modulus and argument of each component of x. By convention we set arg(0) = 0. When we write |x| + c or arg(x) + c, for c ∈ C, the operation is understood componentwise, which means |x j | + c and arg(x j ) + c for all j. Given two vectors x and y in C n , their Euclidean scalar product is denoted by x, y = x * y = x y and means that x, y = n j=1 xj y j , where xj is a complex conjugate. The associated norm is x = x, x 1/2 . Let us consider a matrix A ∈ C m×n , then a j means the j-th row and a •j denotes the j-th column. The component by component product of two vectors is denoted by x y = (x 1 y 1 , . . . , x n y n ) for x, y ∈ C n . The same sign is used if we multiply vector by matrix componentwise which means that for A = (a •1 , . . . , a •n ) ∈ C m×n and x ∈ C m we have that A x = (a •1 x, . . . , a •n x). Let us denote by A(:) the operation of flattening a matrix into a vector such that A(:) = (a 11 , . . . , a 1n , a 21 , . . . , a 2n , . . . , a m1 , . . . , a mn ). Let us also introduce the notations for physical values that are used in this work frequently. By x ∈ C n we denote a set of beams the phases of which we want to correct. The target phases ϕ t = (ϕ t 1 , . . . , ϕ t n ) are encoded into a complex vector x ∈ C n such that x = |x| e iϕ t , where t in superscript means "target". By b ∈ R m + we denote a set of intensity measurements of lasers that were sent through a diffuser. By x ∈ C n we denote an approximation of x by means of some numerical algorithm that uses intensities b for this goal.

Phase retrieval problem

This work studies the phase retrieval algorithms, that can be split into two families:

• The model-based family of algorithms uses a mathematical description of the physical process of diffusion, scattering, and measurements for an array of laser beams. A mathematical model together with measured intensities is then used for phase recovery.

• The model-less family of algorithms tries to approximate the physical transformation by means of a neural network (parametric function) from some family, where a finite set of measurements and phases is required to learn the dependence from scratch.

The conceptual difference between these two families is that the first one uses a model with inputs to produce an output vector, whereas the second case is about building a model that approximates the desired mapping and then can be used for phase recovery. Pictorially, the difference is described on Figure 1.6. These two families are considered in current work as the options for a phase retrieval algorithm block depicted on Figure 1.7 to compare and select the best one.

Let us consider the first family of phase retrieval algorithms. The main idea is to model a signal scattering by a linear map x → { āj , x } m j=1 , where A = (a 1 , . . . , a m ) ∈ C m×n is called a transfer matrix or transmission matrix, n is a number of signals, and m is a number of detectors, and to model signal measurements by a square of modulus x → |x| 2 . Let b 2 ∈ R m + be a vector of intensity measurements returned by physical detectors. Then, the problem of signals reconstruction with a given vector b and known matrix A is called a phase retrieval (see, e.g., [START_REF] Luke | Phase retrieval, what's new?[END_REF][START_REF] Shechtman | Phase retrieval with application to optical imaging: A contemporary overview[END_REF]) and writes as Find x ∈ C n such that | āj , x | = b j , for j ∈ {1, . . . , m},

or equivalently Find x ∈ C n such that |Ax| = b.

(1.2)

1.2. PHASE RETRIEVAL PROBLEM
Note, that a multiplication x by an arbitrary complex number with a unit amplitude does not change the measurements

| āj , x • e iθ | = | āj , x | for all θ ∈ [-π, π[,
which means that a signal x can be reconstructed up to a constant shift in phase.

To use this model a transmission matrix A ∈ C m×n is required. Physically, it depends on the scattering device which is used in the experimental setup. A simple lens corresponds to a Fourier transform of the signal which can be described as the matrix A ∈ C m×n . In this case, knowing the characteristics of a lens the transmission matrix can be easily built. In contrast, if a scattering device is represented by a diffuser, as for the opto-numerical algorithm [START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF], then each j-th row of matrix is computed by means of a phase retrieval algorithm Find a j ∈ C n such that | xk , a j | = b kj for all j ∈ {1, . . . , m}, k ∈ {1, . . . , N },

or equivalently Find A ∈ C m×n such that |XA | = B,
where sets of signals X ∈ C N ×n and measurements B ∈ R N ×m + obtained experimentally are required. To build a transmission matrix A the number of measurements m must be specified in order to have a possibility to solve (1.2). There are results [START_REF] Balan | On signal reconstruction without noisy phase[END_REF][START_REF] Conca | An algebraic characterization of injectivity in phase retrieval[END_REF], where it is proved that when m ≥ 4n -4 then the problem (1.2) always has a unique solution up to a constant phase for any generic A. Following [START_REF] Afonso | Saving phase: Injectivity and stability for phase retrieval[END_REF], the term "generic" can be thought of as the property which is satisfied with probability 1 by measurement vectors {a i } m i=1 ⊆ C n drawn from continuous distributions. In other words, if columns of matrix A were drawn from a continuous probability distribution then mapping x → |Ax| is injective and thus a solution of (1.2) is unique up to a constant phase. Despite the fact of uniqueness, it is still a difficult problem to retrieve vector x ∈ C n by its measurements b ∈ R m + . For this aim, a huge variety of numerical algorithms were developed. The main idea is to reformulate (1.2) in terms of a minimization problem with some useful properties for a specific algorithm. For instance, the reformulation min (x,y)∈C n ×C m 1 2 Ax -y 2 s.t. |y| = b.

(1.3) can be solved in terms of iterative projections between the range of a transmission matrix range

(A) = {y ∈ C m : y = Ax, ∀x ∈ C n } and the set M b = {y ∈ C m : |y| = b}.
This algorithm is known as alternating minimization or alternating projection algorithm [START_REF] Netrapalli | Phase retrieval using alternating minimization[END_REF][START_REF] Marchesini | Alternating projection, ptychographic imaging and phase synchronization[END_REF]. It is also sometimes referred as the Gerchberg-Saxton algorithm [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF], while the original Gerchberg-Saxton algorithm [START_REF] Gerchberg | Practical algorithm for determination of phase from image and diffraction plane pictures[END_REF] is when the matrix A corresponds to a discrete Fourier transform. This formulation can be modified by introducing variables u and θ which replaces y in (1.3) by u e iθ and writes as min

(x,u,θ)∈C n ×R m + ×R m 1 2 u -b 2 s.t. Ax = u e iθ . (1.4)
This reformulation is used in [START_REF] Liang | Phase retrieval via the alternating direction method of multipliers[END_REF] to solve (1.2) by means of Alternating Directions Method of Multipliers [28, §4.4] which is also a part of projection algorithms family. In contrast to the alternating projections algorithm, ADMM requires a way of computing a regularization parameter ρ, which is an additional difficulty.

PHASE CORRECTION PROBLEM

There is also the reformulation which is useful for the descent direction family of algorithms.

min x∈C n 1 2 |Ax| 2 -b 2 2 (1.5)
Using a square of left and right parts of (1.2), the function to minimize in (1.5) is differentiable in terms of Wirtinger calculus (Section 2.3). This property gives an ability to apply the descent algorithms like gradient descent [START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF], Gauss-Newton [START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF] or other.

There is the family of algorithms that is based on the reformulation of (1.2) in terms of a relaxed semi-definite programming [START_REF] Candès | Phase retrieval via matrix completion[END_REF][START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF]. One possible formulation min X∈C n×n trace(X)

s.t. trace(a j a * j X) = b 2 X 0 (1.6)
where X = xx * when exact recovery occurs. Otherwise, it is reasonable to compute a normalized leading eigenvector of X to obtain a solution of (1.2). This relaxation is called as PhaseLift [START_REF] Candès | Phase retrieval via matrix completion[END_REF].

Another reformulation from this family writes as

min U ∈C n×n trace(U M ) s.t. diag(U ) = 1 n U 0 (1.7)
where

M = diag(b)(I -AA † )diag(b). If rank(U ) = 1, then the solution of (1.2) is such that U = xx * .
Otherwise, the normalized leading eigenvector of U is used as an approximate solution of (1.2). This approach known as a MaxCut [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF]. The number of parameters to optimize in (1.7) and (1.6) raise with a quadratic rate if n increases. Since an efficiency of a phase retrieval algorithm is the main criteria in this work, reformulations (1.6) and (1.7) are not considered to use.

Let us consider the second family of phase retrieval algorithms. In this case, there are no assumptions on the form of connection between phases of beams and the measurements. The only information that could be given is a set of

N couples D = {(x i , b i )} N i=1 where x i ∈ C n and b i ∈ R m
+ are input signals and experimentally performed measurements. The goal is to use D as training data to obtain an artificial neural network (NN) that approximates measurements to signals mapping. One of the schemes covered in the published literature [START_REF] Hou | Deep-learning-based phase control method for tiled aperture coherent beam combining systems[END_REF] relies on a phase recovery by a convolutional neural network, such as in the pioneering work on NN for adaptive optics [START_REF] Sandler | Use of a neural network to control an adaptive optics system for an astronomical telescope[END_REF]. The NN serves to map the intensity of an interference pattern of the beam array directly into the distribution of phase in the array. The simulations reported in [START_REF] Hou | Deep-learning-based phase control method for tiled aperture coherent beam combining systems[END_REF] shows that the accuracy of the CNN-based phase control drops when the array increases from 7 to 19 beams. This is a limitation that was also highlighted in the field of wavefront sensing so that NNs were often used only as a preliminary step for initialization of an optimization routine [START_REF] Paine | Machine learning for improved image-based wavefront sensing[END_REF]. The numerical experiments in [START_REF] Wang | Stabilization of the 81-channel coherent beam combination using machine learning[END_REF] also confirms that a direct phase recovery appears to be a complicated task.

Phase correction problem

Schematically, the process of phase control for a set of n beams is described on Figure 1.7. There, one can see that the input beam array x ∈ C n , the phases of which are unknown but with known amplitudes, is passed through the phase modulation block where phases can be changed using a vector ϕ ∈ [-π, π[ n . Then, beams are split by a beamsplitter where a reflected part is sent to the output and a passed part is sent to a diffuser. This device mixes the input beams which produces an interference pattern which is called speckle. The intensities of this pattern are measured by means of photodetectors visualized by green dots on the speckle image and stored in the vector b ∈ R m + , where m > n. These measurements together with the target signal x are then given to the phase retrieval algorithm, the goal of which is to reconstruct approximately the measured signals x using b and x. The phases of a retrieved approximation x ∈ C n are then subtracted from the given phases of the target signal x, and sent to the phase modulation block. The meaning of phase correction vector ϕ = arg( x) -arg( x) can be clearly Figure 1.7: Opto-numerical algorithm [START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF]. seen in the case of one laser (Figure 1.8). The vector ϕ, which is applied to x, moves it closer to the target x. Thus, we have that

arg(x ) = arg(x) + arg( x) -arg( x) = arg(x) + ϕ,
from where it is clear that arg(x ) = arg( x) if retrieved phases arg( x) coincide with arg(x). However, one correction can be not enough to lock the phases to the prescribed values. That is why the process of measurements and corrections repeats several times. At this level, it is important to understand a difference between a phase correction and phase retrieval problems. As it was revealed above, by phase correction task we consider a problem to set beams' phases to the prescribed target values, where the initial beams phases are not observable. A phase retrieval problem aims to approximate x using its measurements b (see, e.g., [START_REF] Luke | Phase retrieval, what's new?[END_REF][START_REF] Shechtman | Phase retrieval with application to optical imaging: A contemporary overview[END_REF]). A particular application of phase retrieval algorithm in a phase correction loop requires target x to be used as an initial point, but originally it is not required.

Noise and data generation

As it is discussed in Chapter 1, the experimental data contains a high level of noise. Thus, one of the main requirements for the algorithms that are considered in this work is robustness. To verify it numerically, the noise must be generated and added during numerical simulations, which gives the ability to estimate the robustness of the algorithm before using it in experiments. This section explains how noise is generated during numerical simulations.

In practice, the noise comes from intensity measurements by means of a camera or photodetectors. Experimental intensity measurements b ∈ R m + together with input signals x ∈ C n are used to find parameters a j ∈ C n for j ∈ {1, . . . , m} of the mathematical model (1.2), which leads to a presence of noise in it. To simulate that numerically it is proposed to use the following approach

A = A + σε, (2.1)
where σ controls the level of noise, ε ∈ C m×n be a random matrix where Re{ε ij }, Im{ε ij } ∼ N (0, 1), A ∈ C m×n be the generated transmission matrix in the same way as ε, which is used to compute a vector of measurements b ∈ R m + , and A is a noisy transmission matrix, which will be given to the algorithm. Thus the algorithm will try to find x ∈ C n such that |Ax| = b using A and b. The same conditions are for initialization methods.

Metrics

In this section, we review several metrics than can be used to evaluate the distance between two signals represented by two vectors x and y in C n .

The first metric, used in [START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF] 1 , is useful to only compare the arguments of two vectors. It is defined by

q(x, y) = 1 - x, y |x|, |y| 2 .
(2.2)

Proposition 1 For all vectors x and y in C n , with all nonzero components, q(x, y) = 0 if and only if there exists c ∈ R such that arg(x) = arg(y) + c.

METRICS

Proof. Let (x, y) ∈ C n × C n such that x k = 0 and y k = 0 for all k = 1, . . . , n. The assertion follows from the triangle inequality 2 . Indeed, we have

| x, y | = n k=1 |x k ||y k |e i(arg(y k )-arg(x k )) ≤ n k=1 |x k ||y k | = |x|, |y| ,
with an equality if and only if arg(y k ) = arg(x k ) + c, with c = arg(y 1 ) -arg(x 1 ). 2

The metric (2.2) is sensitive to changes in |x| and |y| (Example 1) which is redundant for our applications where just phase changes are interesting to consider.

Example 1 (q(x, y) sensitivity) Let us consider x = (1, α•e iϕ ) and y = (1, 1) for some α > 0 and ϕ ∈ [0, 2π]. Then we obtain that

q(x, y) = 1 - x, y |x|, |y| 2 = 1 - |1 + α • e iϕ | 2 |1 + α| 2 = 1 - (1 + α cos ϕ) 2 + (-α sin ϕ) 2 (1 + α) 2 = 1 - 1 + 2α cos ϕ + α 2 (1 + α) 2
Let us visualize this function with respect to α ∈ (0, 10] and ϕ ∈ [0, 2π]. It can be seen on It is proposed to use a normalized quality metric q norm (x, y) = 1 -n -2 | e i arg(x) , e i arg(y) | 2 .

(2.3)

Proposition 2 For all vectors x and y in C n , with all nonzero components, q norm (x, y) = 0 if and only if there exists c ∈ R such that arg(x) = arg(y) + c. Moreover, q norm (x, y) is such that

q norm (x, y) = q x |x| , y |y| 2 .
Proof. The first assertion is shown following the same approach as in Proposition 1 setting |x| = |y| = 1 n . The second assertion is shown by direct formula evaluation for x and y from

C n q x |x| , y |y| 2 = 1 - x |x| , y |y| x |x| , y |y| 2 = 1 - e i arg(x) , e i arg(y) 1 n , 1 n 2 = 1 -n -2 | e i arg(x) , e i arg(y) | 2 2
In [START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF], the distance between two complex vectors is defined as follows. For (x, y)

∈ C n × C n define dist(x, y) = min φ∈[0,2π]
x -e iφ y .

(2.4)

Proposition 3 For all vectors x and y in C n , dist(x, y) = x -e -i arg( x,y ) y . Moreover, dist(x, y) = 0 if and only if |x| = |y| and there exists c ∈ R such that arg(x) = arg(y) + c.

Proof. Let (x, y), ∈ C n × C n and φ ∈ R. We have

x -e iφ y 2 = x 2 + y 2 -2 Re(e iφ x, y ).

The result follows from the fact that Re(e iφ x, y ) ≤ |e iφ x, y | = | x, y |, with an equality if and ony if φ = -arg( x, y ). 2

In some cases, it is also useful to use normalized distance defined as

dist norm (x, y) = dist(x, y) max{ x , y } . (2.5)

Wirtinger calculus

In this section, we explain Wirtinger calculus which is a way to compute derivatives of a realvalued function on a complex-valued domain f : C → R. The application of Wirtinger calculus is used in this work to differentiate metrics that are defined in Section 2.2. It is known from the analysis that complex variable function is differentiable if Cauchy-Riemann conditions are satisfied i.e. the function is holomorphic. However, the neural network model and loss can be non-holomorphic functions. In this case, Wirtinger calculus [START_REF] Kreutz-Delgado | The complex gradient operator and the cr-calculus[END_REF] can be applied. The idea of Wirtinger was that any differentiable mapping f : R 2 → R 2 can be cast in the complex domain by introducing the conjugate coordinates and changing variables such that

z z = 1 i 1 -i x y . (2.6)
It means that any complex-or real-valued function

f (z) = f (x, y) = u(x, y) + iv(x, y),
of several variables can be written in the form f (z, z), where f is holomorphic in z = x + iy for fixed z,and holomorphic in z = x -iy for fixed z. This holds as long as a the real-valued functions u and v are differentiable as functions of the real variables x and y [START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF]. As an example, consider

f (z) = (b -|a * z| 2 ) 2 = (b -z aa * z) 2 = f (z, z), with z, a ∈ C n , b ∈ R. While f (z) is not holomorphic in z, f (z, z
) is holomorphic in z with fixed z and vice versa. Then, the Wirtinger derivative can be defined as in Observation 1.

Observation 1 (Wirtinger derivative) Let us consider a real-valued function f (x, y) : R 2 → R. Let f (z, z) : C 2 → R be the same function defined on a complex domain by changing variables as it is shown in (2.6) for z ∈ C and its conjugate z, and f (x, y) is differentiable with respect to x and y. Then derivatives on the complex domain can be expressed via derivatives in the real domain and write

∂f ∂z := ∂f (z, z) ∂z z=const = 1 2 ∂f (x, y) ∂x -i ∂f (x, y) ∂y ∂f ∂ z := ∂f (z, z) ∂ z z=const = 1 2
∂f (x, y) ∂x + i ∂f (x, y) ∂y .

Proof. From (2.6), we can explicitly express x and y in terms of z and z as

x = 1 2 (z + z), y = - i 2 (z -z). Let us define h x (z, z) = 1 2 (z + z) and h y (z, z) = -i 2 (z -z). Then, since f (z, z) = f (x, y) = f h x (z, z), h y (z, z) , we obtain that ∂f (z, z) ∂z = ∂f (x, y) ∂h x ∂h x ∂z + ∂f (x, y) ∂h y ∂h y ∂z = 1 2 ∂f (x, y) ∂x -i ∂f (x, y) ∂y .
Derivative with respect to z can be obtained in the same way. 2

The idea, presented in Observation 1, can be generalized to a multivariate function in order to compute a Wirtinger gradient (Observation 2). For this aim, let us write the following useful identities, which are applied to f (z, z) for z ∈ C.

∂f ∂z = ∂ f ∂ z , when f is real ∂f ∂z = ∂f ∂ z (Conjugation rule) ∂f ∂ z = ∂ f ∂z , when f is real ∂f ∂ z = ∂f
c = z z = I iI I -iI x y = M r, (2.7)
where I here is an identity matrix n by n, x, y ∈ R n . Note here that M -1 = 1 2 M * . Since the coordinate transformation is a linear transformation it follows that for f :

C n → R that can be represented in different coordinates f (z) = f (z, z) = f (c) = f (r) = f (x, y) we have that ∂f ∂r = ∂c ∂r ∂f ∂c = M ∂f ∂c .
Note, that since f is real-valued then we have that ∂f ∂r = ∂f ∂r = M * ∂f ∂c , which will be used to define gradient in c. Let us define the displacement in r as ∆r = -α ∂f ∂r for some α > 0, which coincides with a displacement in a gradient descent method. Thus, we are interested in how it writes in terms of complex coordinates c. Thus we obtain that

∆c = M ∆r = -αM ∂f ∂r = -αM M * ∂f ∂c = -2α ∂f ∂c (2.8)

TENSORFLOW LIBRARY EXTENSION

Because c = (z * , z ) , (2.8) can be written as

∆c = ∆z ∆z = -2α     ∂f ∂ z ∂f ∂z     .
(2.9)

The proof follows.

2

Note, that from Observation 2 we have that the gradient of function of n complex variables is of dimension 2n. Clearly, it is because the augmented system of coordinates (z, z) is used. Thus, a displacement in the variable z only writes as -2α ∂f ∂ z and then gradient is

∇ z f (z) = 2 ∂f ∂ z .
(2.10)

Tensorflow library extension

Wirtinger calculus can also be applied to compute gradients of complex-valued neural networks. This possibility is already implemented [START_REF] Barrachina | Complex-valued vs. real-valued neural networks for classification perspectives: An example on non-circular data[END_REF] in Tensorflow library, which is a popular framework to build, train, and use neural networks. The gradient of a complex-valued function here is computed as

∇ z f (z) = ∂f ∂z + ∂ f ∂z , (2.11)
which coincides with (2.10) if f is a real-valued function. Despite the fact that Tensorflow supports complex derivatives, there is no support for complex-valued layers. There exist several implementations on GitHub, which extend Tensorflow with complex-valued layers. However, not all of them are maintained by their creators, which produces difficulties in the version control of dependent libraries. For this goal, we created our own implementation (https://gitlab.xlim.fr/shpakovych/cvnn) that extensively uses a polymorphic paradigm in object-oriented programming, which gives an ability to reduce significantly the size of the library in comparison with other implementations.

Profiling method

In this section, we formulate the profiling method to be able to compare the performance of the algorithms, which are considered in this work. The main interest is to measure the required time to solve a problem by a selected algorithm. To achieve this goal we must specify both the conditions under which we conclude that the problem is solved and the way to measure time.

Definition 1 (Solved problem)

The problem is considered to be solved if the value of a distance metric between the k-th approximation of a solution and a solution itself is smaller than a selected tolerance for the last 10% percent of iterations, where the total number of iterations was fixed before. for random numbers 1 million times and then taking the average. For this goal, a program in C language was written.

Example 3 (Time of basic operations) Let us count the time of evaluation sin(x) for x ∈ R n . Following the Table 2.1, we obtain 1.4113816 • 10 -8 n seconds.

Definition 4 (Profiling method) The idea is to build a statistical curve that reveals which part of the problems can be solved by a selected time. The same strategy can be applied for the required number of algorithm iterations to achieve a selected tolerance.

Example 4 (Profiling method) For instance, from Figure 2.2 we can conclude that Alternating projection method with random initialization point can solve approximately 50% of all problems in 4 ms and requires maximum 200 iterations to solve 55% of all problems.

Heat map method

In this section, we introduce the method to explore the capabilities of the selected algorithm. In this work, we consider two kinds of problems to solve: phase retrieval and phase correction. The size of these problems is determined by the number of beams n and the number of intensity measurements m. Also, we consider two kinds of algorithms: optimization methods and neural network methods.

In the context of optimization methods, it is reasonable to check how many problems can be solved for different relations of n and m, and how many iterations it requires. Thus, we generate N problems for different n and m, try to solve them by means of an optimization method with a defined maximal number of iterations and required stopping tolerance, and count how many problems were solved, where the meaning of the solved problem is given in Definition 1. Let us consider that M problems were solved for some fixed n and m. Then, we get the M/N fraction which characterizes the solved part of all problems. Also, we compute how many iterations we require on average to solve the problem of size n by m. Then, it is convenient to visualize these results in terms of heat maps as it is presented on Figure 2.3.

In the context of neural network methods, we are more interested in a potential capability to solve problems of size defined by n and m. That is why instead of the fraction of solved problems we plot an average minimum value of the selected metric in a gray scale. The neural network method does not require iterations to solve a phase retrieval problem. That is why, instead of an average number of iterations, we plot an average relative time required to train the model of size defined by n and m (Figure 2 Chapter 3

Machine learning framework

Let X be a space of objects, Y be a space of responses, and let y * : X → Y be an unknown function that exactly maps each object to its response. Let the evaluations y i = y * (x i ) are known on a given subset {x 1 , . . . , x N } ⊂ X of size N . Let D = {(x i , y i )} N i=1 be a set of couples which we call training sample.

The goal is to find w ∈ W(F) such that a(•, w) ≈ y * (•) for

a : X × W(F) → Y, (3.1)
where a is taken from some parametric family of functions a ∈ {a(x, w) :

X × W(F) → Y }, W (
F) is a space of weights over a field F (either R or C). Each parametric family defines the form of a.

The quality of approximation y * by a depends on the parametric family from where a was taken and on the size of training sample N . Since y * is unknown, we can measure the quality just on the training sample or on some subset of it. Usually, people split the training sample on the train and test disjoint subsets

D train = D N 1 i=1 , D test = D N i=N 1
where N 1 = N -k to avoid evaluation of model on the information that was used to find the parameters w * ∈ W(F). The quality approximation on a given subset can be measured by loss or cost function L : Y ×Y → R.

If a loss function is a differentiable or at least piecewise differentiable then the task of finding parameters w * ∈ W(F) can be formulated as the unconstrained optimization problem w * ∈ argmin w∈W(F ) L(a(•, w), D train ).

(

3.2)

There is a big variety of optimization methods that can be used to solve unconstrained problems. However, this set reduces with an increase of a training sample size and a dimension of parametric space W(F). For example, if it is not possible to store a training sample in the memory or the time of computing gradients is too long, then only stochastic methods like stochastic gradient descent [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] can be applied. These methods have two big advantages that are extensively used in machine learning. First, they can use a descent directions computed just on a subset of a training sample to perform an optimization step. Second, the steplength of a descent direction is computed without linesearch methods which increases an optimization speed dramatically.

When the parameters w * are found, we obtain a * (•) = a(•, w * ) function that approximates y * (•) in sense of loss function on the given training set. The test subset of training sample then is used to check the generalizing ability of the approximation a.

Example 5 Let us consider the regression problem. For instance, we need to know the price y of a flat using the information about a number of rooms x 1 and a distance to the transport line x 2 . And suppose that we know these characteristics for N flats. It means that we have training sample of size

N in form D = {([x 1,i , x 2,i ], y i )} N i=1 where X = {[x 1,i , x 2,i ]} N i=1 is the set of objects and Y = {y i } N
i=1 is the set of responses. Let we also have an assumption that there is a linear dependence between x 1 , x 2 and y. It means that we can choose a : X × W(F) → Y from the parametric family of linear functions

a(x, w) = w 0 + w 1 x 1 + w 2 x 2 , where x = [x 1 , x 2 ] ∈ R 2 , w = [w 0 , w 1 , w 2 ] ∈ W(R) = R 3 .
Then we define the loss function for this task. The obvious choice is the squared l 2 norm normalized by the number of objects N

L(ŷ, y) = 1 N N i=1 (ŷ i -y i ) 2 2 ,
where ŷi = w 0 + w 1 x 1,i + w 2 x 2,i is the predicted i-th price, and y i is the original i-th price from training sample. Then, we split training sample on the train and test disjoint subsets D train and D test and define X train , Y train to be objects and responses from D train , and X test and Y test to be objects and responses from D test . Then the problem of finding parameters w can be formulated as the following optimization problem

w * ∈ argmin w∈W(R) 1 |D train | |D train | i=1 a([X train ] i , w) -[Y train ] i 2 .
After finding the optimal parameters we can evaluate this model on the test set

score = 1 |D test | |Dtest| i=1 a([X test ] i , w) -[Y test ] i 2 ,
and finish.

In the example 5, the step-by-step process of building and training the simple linear model with a squared error loss function is presented. The linear parametric family of functions which is used in 5 is appropriate for this simple problem, however its complexity increases with the complexity a problem.

Typically, a parametric family is a composition of linear blocks with nonlinear transformations between them. A common name for the elements that are used for building a parametric family is layer. There are two types of layers: trainable and non-trainable. By trainable we mean that it contains parameters to be optimized during learning. In example 5 it was the vector w = [w 0 , w 1 , w 2 ]. In the opposite case, there are no parameters to optimize i.e. some data transformation like dropout regularization (see Section 3.3).

LAYERS

Layers

In this section, we consider an example of trainable layer which is called fully connected or dense layer. The fully connected layer is an affine function fc : The optimization requires gradients with respect to trainable variables. In the case where F = R we have

F n × F m×n × F m → F m that writes as fc(x, W, b) = W x + b, (3.3
∇ W fc(x, W, b) = [x, x, . . . , x] ∈ R m×n , ∇ b fc(x, W, b) = 1 R m , (3.4)
where 1 R m is a vector of ones of dimension m. In the case where F = C we have

∇ W fc(x, W, b) = [x, x, . . . , x] * ∈ C m×n , ∇ b fc(x, W, b) = 1 C m , (3.5)
where 1 C m is a vector of 1+0i elements of dimension m. The complex differentiation is explained in Section 2.3.

Let us give a simple example of a binary classification problem that can be solved by one layer model. The only thing that we add which was not discussed before is a sigmoid activation function. It scales the output of a fully connected layer into the range from 0 to 1. we give a label 0 for objects x that were sampled from N (-2, 1) and 1 for x that were sampled from N (2, 1).

Thus Y = [Y 1 , Y 2 ]
× R 2 × R → R fc(x, w, b) = w 1 x 1 + w 2 x 2 + b, (3.6)
where x ∈ R 2 is the input, w ∈ R 2 and b ∈ R are the parameters to optimize. Then, since our responses are ether 0 or 1 we will scale the output of fc(x, w, b) into the range [0, 1] using a sigmoid function.

σ(x) = 1 1 + exp(-x) . (3.7)
Thus the parametric family writes as

a(x, w, b) = 1 1 + exp(-(w 1 x 1 + w 2 x 2 + b)) . (3.8)
Then we need to define the loss function. For a task of binary classification where responses are either 0 or 1 there is well-known loss function that is called log loss of binary cross-entropy that came from the information theory and writes as However, log loss is difficult to interpret in terms of quality of the approximation. In this case, we need to compute some additional metrics like accuracy, which is simply a fraction of the correct predictions to all samples in a set. The learning process is visualized in Figure 3.2 with log loss and accuracy metric on the vertical axis and the number of gradient steps on the horizontal axis. After finding the optimal parameters we can evaluate this model on the test set

L(ŷ, y) = - 1 N N i=1 y i log(ŷ i ) + (1 -y i ) log(1 -ŷi ). ( 3

Activation functions

There is a wide range of special functions that can be applied for model building. Usually, they are added between layers that allow learning nonlinear hidden dependencies in data. Even though the number of activation functions exists, some core functions are commonly used. One of then is a sigmoid function that was used in Example 6. It scales the input data into the range from 0 to 1 by the formula

σ(x) = 1 1 + exp(-x) ,
and the plot is depicted on Figure 3.4a. Usually, this layer is applied as a last one and mainly

REGULARIZATION

There are also some modifications of ReLU that can be useful for better convergence. They are LeakyReLU and ELU (exponential linear unit).

LeakyRelu(x, α) = αx x < 0 x x ≥ 0 , ELU(x, α) = α * (exp(x) -1) x < 0 x x ≥ 0 (

Regularization

Regularization is a layer for reducing the overfitting (or memorization) phenomenon during learning. The indicator of overfitting is a different behavior of loss value during the learning process on a train and test sets. The case when loss decreases on a train set and increases on a test is called overfitting. The classical regularizers are l 1 or l 2 norm of the trainable parameters, which is added as an additional term to a loss function to prevents high variance for model parameters that leads to its better generalization. However, for neural networks, there is a more efficient state-of-the-art regularizer which is called dropout [START_REF] Srivastava | Improving neural networks by preventing co-adaptation of feature detectors[END_REF]. The idea is randomly set to zero some parts of neurons during training. This encourages each neuron to be independently useful and does not rely on the output of other neurons. More precisely, with some predefined probability p ∈ [0, 1] the elements of vector y (which is for example y = fc(x)) are set to zero. Other elements are scaled up by coefficient 1 1-p so that the expected value is preserved.

Loss function

The loss function or cost function is a function that maps input values I ⊂ R n onto a real number which can be interpreted as a cost.

L : I → R. (3.12)
This kind of functions is used to measure a quality of approximation of an unknown dependence between inputs and outputs by some model a. Each kind of tasks requires a special loss function.

For instance, mean squared error measure is useful for regression task

3.5. GRADIENTS COMPUTING MSE(ŷ, y) = ŷ -y 2 2 , (3.13)
where ŷ ∈ R n is predicted values by a model, y ∈ R n is the original values that were taken from the training sample. However, this function starts being less useful for classification tasks for m classes where a cross-entropy loss plays an important role and writes as

H(ŷ, y) = - 1 n n k=1 m i=1 y k,i log(ŷ k,i ), (3.14)
where y, ŷ ∈ [0, 1] n×m are predicted and true classes where y k ∈ R m is a one-hot vector (full of zeros except one entry with 1) and ŷk ∈ R m is a probability vector (ŷ k ≥ 0 for k ∈ {1, . . . , m} and m i=1 ŷk,i = 1).

Gradients computing

The step of computing gradient is crucial for understanding the learning process. As was mentioned at the beginning of this chapter, by learning we consider a process of minimization of a loss function with respect to trainable parameters. In most of the cases, machine learning uses the first-order derivatives for optimization. The main motivation is a computational speed which is a crucial part in machine learning. Usually, the model is represented by the composition of the functions which leads us to the result in analysis about differentiation of a composite function that is known as a chain rule.

d(f 1 • • • • • f n )(x) dx = df 1 df 2 • • • • • df n dx .
This result is a core for model differentiation in a machine learning world which is used to define a Back-propagation algorithm, which is just the another name for algorithmic differentiation [START_REF] Günş Baydin | Automatic differentiation in machine learning: a survey[END_REF] method that was developed to track an inner data interaction and efficiently compute derivatives. We will not give a detailed explanation of this area in current work but just consider that gradients at given point x 0 ∈ F n for a composite functions can be computed fast for millions of trainable parameters.

Supervised learning

When the training sample is collected, the model is built, the loss function is defined we are ready to train or learn the model. Recall, that by learning, we mean the process of solving an optimization problem of the form

w * ∈ argmin w∈W(F ) 1 N N j=1 L(a(x j , w), y j ), (3.15)
where W(F) is a space of trainable parameters over field F (R or C), N is the size of training sample {(x i , y i )} N i=1 , x i and y i for i ∈ {1, . . . , N } are objects and responses, a(x, w) is the function from some parametric family. Then, learning algorithm 1 can be applied. (e) Let U : W(F) × W(F) → W(F) is a functions that computes parameters update with respect to a given descent direction ∆w and state w.

Set k = 0.

2. Gradients computing:

∆w k = ∇ w L(a(X train , w), Y train )| w=w k 3. Gradients applying: w k+1 = U(w k , ∆w k ) 4.
Return loop: Set k = k + 1 and go to 2.

Quasi-Reinforcement learning

In this section, we consider a special kind of learning, which is called Reinforcement learning or simply RL. To be more precise, we need discuss a modification of this process that has been developed in this work. However, to have a full picture in mind, we need to start from the original algorithm. High-level definition reveals that the reinforcement learning is devoted to learn an agent how to do an action a t ∈ A at the current state s t ∈ S to maximize reward r t ∈ R of the environment E : A → S × R . Let us give a simple example of the reinforcement learning problem where all of these parts are presented.

Example 7 Let as consider a stick balancing problem. A stick is attached by an un-actuated joint to a cart, which moves along a frictionless track. This is our environment E : A → S × R. The characteristics of the stick at time t are 4 numbers: linear and angular position and velocity. They are the state and thus s t ∈ S = R 4 . The system is controlled by applying a force of +1 or -1 to the cart. The ±1 force is the action a t ∈ A = {-1, 1}. A reward r t ∈ R = {+1} is provided for every timestep that the stick remains upright. The goal is to create agent : R 4 → {-1, 1} that keep the stick from falling over and thus maximize the total reward r = T t=0 r t where T is a maximal number of steps per simulation.

Reinforcement learning approach proposes to find the agent from the parametric family of neural network functions. The main difficulty is that we can not formulate it as supervised problem since we do not know the correct action a t at state s t and thus the reward is not differentiable because it computes by environment which is considered as a black box. This fact leads to the different approaches of learning agent which are all based on the repetitive collection of a huge number of triplets {(a t , s t , r t )} T t=0 that are used to update weights of an agent. However, in this section, we are not concentrated on the explanation of the ways of updating weights using the triplets, but on the problem when correct actions are known and can be returned from the environment E : A × A → S × A, and, as a consequence, the reward can be calculated in the same way as a loss. The motivation is that the phase correction problem can be formulated in this framework. Thus we call it quasi-reinforcement learning to point on the fact that it is not classical case. The algorithm of quasi-reinforcement learning can be formulated as the Algorithm 2. Note, that it requires a possibility to easily generate states an correct actions at step 2 to perform learning. The application of this approach to the phase correction problem is presented in Chapter 6.

Algorithm 2: QRL algorithm 1. Initialization:

(a) Let w (0) ∈ W(F) be the starting value for trainable parameters.

(b) Let Agent : S × W(F) → A be a model to learn.

(c) Let R : A × A → R is a reward function. (d) Let E : A × A → S × A is an environment function.
(e) Let U : W(F) × W(F) → W(F) computes parameters update concerning a given descent direction ∆w and state w.

(f) T ∈ N\{0} is the maximal number of steps and N ∈ N\{0} is the size of training sample.

Set k = 0.

Training samples generation:

To generate initial states and correct actions

S 0 = {s (i) 0 } N i=0 , A 0 = {a (i) 0 } N i=0 .
3. Environment discovering: Set t = 0, w

(k) t = w (k) .
(a) To predict actions Ât = Agent(S t , w

(k) t ). (b) To compute gradients ∆w (k) t = ∇ w R( Ât , A t )| w=w (k) t (c) To apply gradients w (k) t+1 = U(w (k) t , ∆w (k) t ) (d) To send actions to the environment S t+1 , A t+1 = E( Ât , A t ) (e) Set t = t + 1. If t < T then w (k) = w (k) t
and go to (3a).

Chapter 4

Phase retrieval algorithms

In this chapter we present a set of phase retrieval algorithms that can be used as a numerical part of the opto-numerical algorithm [START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF]. The first two sections Section 4.2 and Section 4.3 describe the descent direction optimization algorithms: gradient descent and Gauss-Newton. It is known that the gradient descent iteration is cheep if a line search strategy selected properly. It makes this algorithm a good candidate for the application in a phase correction loop. Also, despite the fact that the Gauss-Newton method requires just the first order derivatives, the calculation of a Jacobian and solving the system is a computationally expensive part. However, it is known that this algorithm requires far less iterations to achieve a selected tolerance. This is the reason to check this algorithm also. Then, Section 4.5 and Section 4.6 describe the optimization algorithms from the family of the projection methods: the Alternating Projections (AP) and the Alternating Direction Method of Multipliers (ADMM). It is known that AP and ADMM require a low number of operations for one update, which also makes them good candidates for current applications. Another point that is considered in this chapter and discovered in Section 4.7 is the initialization strategies [START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF][START_REF] Liu | Phase retrieval via wirtinger flow algorithm and its variants[END_REF] which aim to find a starting point in the neighborhood of a solution. It is shown in this work that a good initial point can increase a speed of the phase correction algorithm significantly. The last section presents the numerical results which show the capabilities of the considered algorithms.

Research path

The goal that we tried to achieve in this section is to try different optimization algorithms to solve a phase retrieval problem and to find the best one with respect to the requirements defined in the introduction section. As was mentioned before, there are many reformulations of the original problem which leads to different phase retrieval algorithms. So the goal was also to find a proper formulation.

The most interesting result we achieved was with the ADMM optimization algorithm. We started from the publication [START_REF] Liang | Phase retrieval via the alternating direction method of multipliers[END_REF] where the authors implement the ADMM algorithm to solve a phase retrieval problem. The crucial thing that was observed is that if we decrease ρ we increase the number of problems that the algorithm can solve. However, there was no possibility to set ρ = 0 because of the division operation. That is why we tried to adapt update formulas in order to remove ρ completely. In this way, we discovered the Algorithm 8. Then, we searched for the minimization problem which corresponds to these update formulas and found (4.21). The next step was to ensure that the algorithm is robust. Unfortunately, even with a small level of noise the algorithm does not converge (see Figure 4.8b). This observation can be easily explained since the system (1.2) is overdetermined, which implies that, almost surely, there does not exist y ∈ range(A) and z ∈ M b such that y = z. As a result there is no limit point for the sequence λ (k+1) = λ (k) + y (k+1) -z (k+1) (see Proposition 8).

To overcome this difficulty we introduced the relaxed minimization problem (4.22), which aims to replace the strong linear constraint by y -z = ξ. This gives us the optimization algorithm Algorithm 9. However, then we need to deal with a regularization parameter ρ. The formulas for the adaptive update of ρ can be found in (4.25) and (4.26). Here, the idea is simple. When we are far from a solution, ρ must be small to guide the sequences z (k) and y (k) to the "perfect solution" even if it does not exist. As soon as the distance between z (k) and y (k) is small, then we set ρ = 1, which means that we are using the alternating projection algorithm.

At the end of Section 4.6, we found the algorithm that satisfies all the requirements and, at the same time, converges faster than all other algorithms considered in this chapter.

Gradient Descent Method

In this section we describe the gradient descent method for solving the optimization problem

min x∈C n f (x) := 1 2m |Ax| 2 -b 2 2 . ( 4.1)
It can be easily seen that the solution of (4.1) coincide with the solution of (1.2). Using the Wirtinger calculus (2.10) it is possible to compute the gradient, which in a matrix form writes as

∇ x f (x) = 1 m A * (|Ax| 2 -b 2 ) Ax ,
and can be used in the gradient descent step x (k+1) = x (k) -α∇ x f (x (k) ) where α ∈ (0, 1] is a step length to be computed by one of the line search algorithms and x (0) ∈ C n is an initial point for the algorithm. Then, the gradient descent method can be formulated as Algorithm 3. Algorithm 3: Gradient descent algorithm.

Input: Initial point x (0) ∈ C n Output: Approximate solution x ∈ C n of (1.2). 1 Set k = 0; 2 Loop 3 p (k) = -∇ x f (x (k) ); 4 compute a step length α (k) by some line search method; 5 x (k+1) = x (k) + α (k) p (k) ; 6 if stopping test is satisfied then 7 break 8 k = k + 1; 9 return: x = x (k+1)
The most computationally expensive part is the process of finding the step length α (k) . Thus, it is extremely important to select a line search strategy appropriately in order to maintain efficiency of the Algorithm 3. One of the simplest line search approach is a backtracking algorithm [START_REF] Nocedal | Numerical Optimization[END_REF], which aims to find α that satisfies the following condition

f (x) > f (x -α∇ x f (x)) + c 1 α ∇ x f (x) 2 , (4.2)
where c 1 ∈ (0, 1). The idea is simple: starting from α = 1, we decrease it with some rate until (4.2) is satisfied. This approach is presented in Algorithm 4.

Algorithm 4: Back-tracking line search [START_REF] Nocedal | Numerical Optimization[END_REF].

Input: Initial step length α (0) ∈ (0, 1), current point x ∈ C n , gradient ∇ x f (x) evaluated in x, decrease rate ρ ∈ (0, 1), constant c 1 ∈ (0, 1) Output: Step length α that satisfies (4.2). 1 k = 0; 2 while f (x) < f (x -α (k) ∇ x f (x)) + c 1 α (k) ∇ x f (x) 2 do 3 α (k+1) = ρα (k) ; 4 k = k + 1; 5 return: α = α (k)
Let us provide numerical simulations to observe how Algorithm 3 can solve phase retrieval problem (1.2) for n = 8 and m = 4n with Algorithm 4 in step 4. The results of numerical simulations are presented on Figure 4.1, where the metric dist norm between approximation x (k) and solution is drawn by gray lines for 100 problems. The main interest is to see the behavior of the algorithm with and without the presence of noise. The thing to pay attention is how the value of metric behaves when there is no possibility to improve an approximation due to the presence of noise in the model. For this goal, different solving tolerance is used: 0.001 for Figure 4.1a and 0.2 for Figure 4.1b. The second value was selected empirically from the behavior of the algorithm. It can be seen on Figure 4.1b that there is a limit in dist norm for all of the problems because of the presence of noise in the mathematical model. We observe no oscillations when the metric value reaches its limit, which indicates a sufficient robustness to the noise from this point of view. Also, in the case of no noise, Algorithm 3 with Algorithm 4 in step 4 can solve 95% of all problems and with noise σ = 0.1 it can solve 81% of all problems. Recall, that the problem is considered to be solved with respect to a Definition 1. Thus, there is difference in 14% of solved problems between noisy model and the model with no noise. Also, the starting point was selected randomly but the same for each experiment: with and without noise.

Note that in Algorithm 4 we require evaluation of f at each k-th step in order to check the stopping condition. It involves matrix multiplication at each step, which is extremely expensive for the current application. That is why, let us consider a Two-Point Step Size line search by Jonathan Barzilai and Jonathan Borwein [START_REF] Barzilai | Two-Point Step Size Gradient Methods[END_REF]. This approach is based on the usage of secant equation for hessian approximation of the form αI where I is an identity matrix. According to [START_REF] Barzilai | Two-Point Step Size Gradient Methods[END_REF] the step length α can be computed from the following two optimization problems

min α∈R s (k) -αy (k) 2 , min α∈R αs (k) -y (k) 2 , (4.3)
where

s (k) = x (k) -x (k+1) and y (k) = ∇ x f (x (k) ) -∇ x f (x (k-1)
). The solutions can be easily obtained and written as

α = 2 s (k) , s (k) s (k) , y (k) + y (k) , s (k) , (4.4
α = s (k) , y (k) + y (k) , s (k) 2 y (k) , y (k) , (4.5)
respectively. However, this step length can be applied just when s (k) , y (k) + y (k) , s (k) > 0 which is an analogue to the curvature condition generalized for the case of complex variable. The case of negative curvature is possible since the function to optimize is not convex and if it is detected then backtracking line search is used. Then, we can formalize gradient descent method with the secant based line search in the Algorithm 5. Algorithm 5: Gradient descent algorithm with secant based line search.

Input: Initial point x 0 ∈ C n , decrease rate ρ ∈ (0, 1), constant c 1 ∈ (0, 1). Output: Approximate solution x ∈ C n of (1.2). 1 k = 0; 2 Loop 3 g k = ∇ x f (x k ); 4 if k = 0 or s (k) , y (k) + y (k) , s (k) < 0 then 5 compute α k using Algorithm 4 ;
// use Armojo if negative curvature 6 else 7 compute α k using (4.4) or (4.5) is clear that Algorithm 5 performs faster than Algorithm 3 with Algorithm 4 in step 4. The number of problems that can be solved coincides. Also, we observe no oscillations when it reaches the limit in improvement. To conclude about required time to solve problems by these two algorithm let us build time profiles as it is explained in Section 2.5. It can be seen that in both cases, without noise and with noise, line search based on a secant equation is extremely faster than standard backtracking line search, which is also revealed on

8 x (k+1) = x (k) -α (k) ∇ x f (x (k) );

Gauss-Newton Method

In this section we describe the Gauss-Newton method for solving the optimization problem

min x∈C n f (x) := 1 2m |Ax| 2 -b 2 2 , (4.6)
which is clearly equivalent to (1.2). Following the general scheme of Gauss-Newton method [START_REF] Nocedal | Numerical Optimization[END_REF], we denote a residual function r :

C n → R m + as r(x) = |Ax| 2 -b 2 , (4.7)
and its jacobian matrix in terms of Wirtinger calculus (Observation 2) writes as

∇r(x) = A Āx Ā Ax ⊂ C m×2n . (4.8)
Then, the descent direction p ∈ C 2n is a solution of the system ∇r(x) * ∇r(x)p = -∇r(x) * r(x). (4.9)

However, it can be shown (Observation 3) that the matrix ∇r(x) * ∇r(x) is always not a full rank matrix. Thus, the regularization is required which simply adds δI to ∇r(x) * ∇r(x), where δ > 0. Typical value of δ depends on the data type that is used to represent numbers. For instance, if real and imaginary parts of complex numbers in ∇r(x) * ∇r(x) are represented as float numbers where each takes 32 bits in the memory, then δ must be greater than 10 -8 to provide a regularization effect.

Observation 3 Let m > 2n, then the rank of the jacobian matrix (4.8) is always less than 2n for any x ∈ C n .

Proof. To show that it is enough to find a nonzero vector v ∈ C n such that ∇r(x)v = 0 for any

x ∈ C n . Let v = (x , -x ) then A Āx Ā Ax x -x = (A Āx)x -( Ā Ax)x = Ax Āx -Āx Ax = |Ax| 2 -|Ax| 2 = 0.
The proof follows. 2

Then, system (4.9) is replaced by

∇r(x) * ∇r(x) + δI p = -∇r(x) * r(x), (4.10)
where δ > 0 and I is an identity matrix of size 2n by 2n. From Figure 4.4 and Figure 4.5 we can conclude that both line search strategies are equivalent since the same behavior and the same number of solved problems (≈ 80%) are observed. In addition, no oscillations are found when the algorithm reaches a limit of approximation in the noisy case. Also, from Figure 4.6 we can conclude that both algorithms perform with the same computational speed.

Semidefinite Programming Method

In this section, we consider the convex relaxation of (1.2), which is called PhaseLift and reformulates the original system of nonlinear equations as a semidefinite programming problem (Candes et al. [START_REF] Candès | Phase retrieval via matrix completion[END_REF]) The reformulation is based on a fact that the intensity measurements in (1.2) can be written as

|a * j x| 2 = tr(x * a j a * j x) = tr(a j a * j xx * ) = tr(A j X), (4.11)
where A j is a rank-one matrix a j a * j . Let A(X) = {tr(A j X) : j ∈ {1, . . . , m}}, then the problem (1.2) can be written as find

X s.t. A(X) = b X 0 rank(X) = 1, (4.12)
or equivalently min

X∈C n×n rank(X) s.t. A(X) = b X 0, (4.13)
which appears to be a rank minimization problem. The problem (4.13) is NP hard, and thus it was proposed in [START_REF] Candès | Phase retrieval via matrix completion[END_REF] to use the trace norm as a convex surrogate [START_REF] Beck | Computational study and comparisons of lft reducibility methods[END_REF][START_REF] Mesbahi | On the rank minimization problem over a positive semidefinite linear matrix inequality[END_REF] for the rank functional, which gives the following SDP problem min X∈C n×n tr(X)

s.t. A(X) = b X 0. (4.14)
To recover an original solution x the leading eigenvector of X can be used. Recall, that one of the main requirements to the phase retrieval algorithm defined in Chapter 1 is an efficiency. The number of parameters in reformulation (4.14) is n 2 , where n is the number of laser beams. Since in practice, the number of beams can be up to 100, this method will not be considered further in this work due to obvious inconsistency with the efficiency requirement.

Alternating Projections Method

In this section, we present the numerical solution of (1.2) by means of the alternating projection algorithm. Let us first reformulate equation (1.2) under the equivalent form

Ax = y |y| = b, (4.15)
where the unknown is the couple (x, y) ∈ C n × C m . Obviously, x is a solution of (1.2) if and only if (x, b e i arg(Ax) ) is a solution of (4.15). The system of equations (4.15) can be reformulated as a minimization of a quadratic convex function subject to a nonconvex set of equality constraints:

min (x,y)∈C n ×C m 1 2 Ax -y 2 s.t. |y| = b. (4.16)
We have the following elementary properties.

Proposition 4 Let x be fixed vector in (4.16), then the optimization problem is separable with respect to y j for j ∈ {1, . . . , m}. Using the fact that the constrains are defined for each variable y j separately we can conclude that the minimization of a sum in (4.17) is equivalent to m minimization problems of the form

min y j ∈C |α j -y j | 2 s.t. |y j | = b j , (4.18)
for j ∈ {1, . . . , m} 2

Proposition 5 x ∈ C n is a solution of (1.2) if and only if (x, b e i arg(Ax) ) is an optimal solution of (4.16). In addition, when y ∈ C m is fixed, (4.16) is a linear least squares problem, those optimal solution is x = A † y. When x ∈ C n is fixed, the optimal solution of (4.16) is the orthogonal projection of Ax on the set of constraints and is given by y = b e i arg(Ax) .

Proof. The first assertion follows from the equivalence between (1.2) and (4.15). When y is fixed in (4.16), the constraints disappears and thus the problem is reduced to a linear least squares problem.

When x is fixed in (4.16), following the Proposition 4, the problem is separable and is reduced to the solution of m one dimensional optimization problems of the form (4.18) From the triangle inequality, we know that for all θ ∈ R, |α j -b j e iθ | ≥ | |α j | -b j |, with an equality if and only if θ = arg(α j ) or α j = 0 or b j = 0. The optimal y follows.

2

Proposition 6 For all (x, y) ∈ C n × C n , we have |x| -|y| ≤ x -y , with an equality if and only if for all k ∈ {1, . . . , n}, x (k) = 0 or y (k) = 0 or arg(x (k) ) = arg(y (k) ).

Proof. Let (x, y) ∈ C n × C n . Using the fact that x = |x| , we have

|x| -|y| 2 = x 2 -2 |x|, |y| + y 2 and x -y 2 = x 2 -2 Re( x, y ) + y 2
We have Re( x, y

) = n k=1 |x (k) ||y (k) | cos(arg(y (k) ) -arg(x (k) )) Re( x, y ) = n k=1 |x (k) ||y (k) | cos(arg(y (k) ) -arg(x (k) )) ≤ n k=1 |x (k) ||y (k) | = |x|, |y| ,
with an equality if and only if

n k=1 |x (k) ||y (k) |(1 -cos(arg(y (k) ) -arg(x (k) ))) = 0.
Because each element of the sum is nonnegative, this occurs if and only only one of the modulus or the cosine is zero for all k.

2

Unfortunately, the problem (4.16) is not convex. To handle this difficulty, a well-known approach, suggested by Proposition 5, is to alternatively minimize in x and in y. This leads to Algorithm 7. This algorithm is known as alternating minimization or alternating projection algorithm (see, e.g., [START_REF] Netrapalli | Phase retrieval using alternating minimization[END_REF][START_REF] Marchesini | Alternating projection, ptychographic imaging and phase synchronization[END_REF]). It is also sometimes referred as the Gerchberg-Saxton algorithm (see, e.g., [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF]), while the original Gerchberg-Saxton algorithm [START_REF] Gerchberg | Practical algorithm for determination of phase from image and diffraction plane pictures[END_REF] is when the matrix A corresponds to a discrete Fourier transform (DFT) and the couple of iterates is defined by means of the recurrence formulas y (k+1) = b e i arg(Ax (k) ) and x (k+1) = a e i arg(A † y (k+1) ) , where a is the known vector of amplitudes in the initial space, see e.g., [START_REF] Shechtman | Phase retrieval with application to optical imaging: A contemporary overview[END_REF]. Algorithm 7 shares a Algorithm 7: Alternating projection algorithm.

Input:

Initial point x (0) ∈ C n Output: Approximate solution x ∈ C n of (1.2). 1 k = 0; 2 Loop 3 y (k+1) = b exp(i arg(Ax (k) )); 4 x (k+1) = A † y (k+1) ; 5 k = k + 1; 6 if stopping test is satisfied then 7 break 8 return: x = x (k) .
common property with the original Gerchberg-Saxton algorithm, which is the reduction of the norm of the residual of (1.2) [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF]. Unfortunately, the norm of the residual in Proposition 7 is nonincreasing, which means that Algorithm 7 is not guaranteed to converge to a global minimum of (4.16). Proposition 7 At each iteration k ≥ 0 of Algorithm 7, we have

Ax (k+1) -y (k+1) ≤ Ax (k) -y (k) and |Ax (k+1) | -b ≤ |Ax (k) | -b .
Proof. Let k ∈ N be the iteration number. By using the fact that x (k+1) is a minimum of problem (4.16) with y = y (k+1) and the fact that y (k+1) is a minimum of problem (4.16) with x = x (k) , we have

Ax (k+1) -y (k+1) ≤ Ax (k) -y (k+1) (4.19) ≤ Ax (k) -y (k) ,
which proves the first inequality. To prove the second inequality, use Proposition 6 and (4.19), then apply again Proposition 6 and the fact that arg(y (k+1) ) = arg(Ax (k) ), so that

|Ax (k+1) | -b = |Ax (k+1) | -|y (k+1) | ≤ Ax (k+1) -y (k+1) ≤ Ax (k) -y (k+1) = |Ax (k) | -|y (k+1) | = |Ax (k) | -b . 2
Let us provide numerical simulations to test the robustness of the Algorithm 7. Recall, that the main interest is to see the behavior of the algorithm with and without the presence of noise. The thing to pay attention is how the value of metric behaves when there is no possibility to improve an approximation due to the presence of noise in the model. The results are presented on can conclude that no oscillations are found when the algorithm reaches a limit of approximation in the noisy case. Also, approximately the same number of problems can be solved in both cases: 78% when there is no noise and 74% with the presence of noise.

Alternating Directions Method of Multipliers

In this section, we present the numerical solution of (1.2) by means of the Alternating Directions Method of Multipliers. Let us first reformulate equation (1.2) under the equivalent form

y = z y ∈ range(A), z ∈ M b , (4.20)
where

M b = {z ∈ C m : |z| = b}
z (k+1) = b exp i arg(y (k) + λ (k) ) , y (k+1) = AA † (z (k+1) -λ (k) ), λ (k+1) = λ (k) + y (k+1) -z (k+1) , where A † = (A * A) -1 A * is a pseudo-inverse of matrix A ∈ C m×n .
Proof. The scaled augmented Lagrangian writes as

L ρ (y, z, λ) = Re{ λ * (y -z)} + ρ 2 y -z 2 = ρ 2 y -z + λ 2 - ρ 2 λ 2 ,
where λ = λ ρ is a scaled dual variable. Then, we conclude that 1. The update formula for z with fixed y = y (k) and λ = λ (k) writes as

z (k+1) = argmin z∈M b ρ 2 y (k) -z + λ (k) 2 = b exp i • argmin θ∈R m ρ 2 y (k) -b e iθ + λ (k) = b exp i • arg(y (k) + λ (k) ) .
2. The update formula for y with fixed z = z (k+1) , λ = λ (k) writes as

y (k+1) = argmin y∈range(A) ρ 2 y -z (k+1) + λ (k) 2 = A argmin x∈C n ρ 2 Ax -z (k+1) + λ (k) 2 = AA † (z (k+1) -λ (k) ).
3. The update formula for λ with fixed z = z (k+1) and y = y (k+1) writes as

λ (k+1) = λ (k) + y (k+1) -z (k+1) .
2

The Observation 4 is used for simplification the update formula for y (k+1) .

Observation 4 Let us consider the update formulas in Proposition 8. If A * λ (0) = 0 then

A * λ (k) = 0 for any k ∈ N\{0}.
Proof. Let us prove the result by induction reasoning. By assumption the result is true for k = 0. Suppose that A * x (k) = 0 for k ∈ N. We have

λ (k+1) = λ (k) + AA † z (k+1) -z (k+1) we have that A * λ (k+1) = A * λ (k) + A * z (k+1) -A * z (k+1) and thus A * λ (k+1) = 0. 2
Algorithm 8: Alternating Directions Method of Multipliers.

Input: Initial point x (0) ∈ C n and parameter λ (0) ∈ C m Output: Approximate solution x ∈ C n of (1.2). 1 k = 0; 2 y 0 = Ax (0) ; 3 Loop 4 z (k+1) = b exp(i arg(y (k) + λ (k) )); 5 x (k+1) = A † z (k+1) ; 6 y (k+1) = Ax (k+1) ; 7 if stopping test is satisfied then 8 break 9 λ (k+1) = λ (k) + y (k+1) -z (k+1) ; 10 k = k + 1; 11 return: x = x (k) .
Proposition 9 Let (y, z, λ) ∈ C m × C m × C m be a limit point of a sequence defined by update formulas in Proposition 8 and A * λ (0) = 0, then (y, z) is a global solution of (4.21).

Proof. Let (y, z, λ) ∈ C m × C m × C m be a limit point of a sequence defined by update formulas in Proposition 8 exists and A * λ (0) = 0. Then, replacing y by AA † z we have that

z = b exp i • arg(AA † z + λ) , λ = λ + (AA † -I)z.
We can immediately deduce that z ∈ M b and (AA † -I)z = 0. The matrix AA † -I is a projection matrix on the kernel of A * . The projection is zero only if z ∈ range(A).

2

Remark 1 By Proposition 9 and Proposition 8 we conclude that if the Alernating Directions Method of Multipliers converges, then it always returns a global solution of (4.21) and so (1.2). This is an important property which the alternating projections algorithm does not have [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF].

In practice, it is possible that there is no exact solution of (1.2) because of noisy measurements b and inaccurate transmission matrix A. In context of problem (4.21), it means that range(A) ∩ M b = ∅. Note, that this leads to unpleasant consequence for update formula of dual multipliers λ in Proposition 8. It is clear that if range(A) ∩ M b = ∅, then there are no y ∈ range(A) and z ∈ M b such that y -z = 0 and thus no limit point for λ.

Let us provide numerical simulations to demonstrate this fact. On Figure 4.8 it can be clearly seen that when there is noise in the model with σ = 0.1, then Algorithm 8 does not converge. However, when no noise is added to the model, then all of 100 problems can be solved by means of Algorithm 8 despite the fact of extensive oscillations during convergence. This fact has been never observed for previous methods discussed in this chapter. Since the presence of 0 100 200 300 400 500 Iterations noise is true for the applications that are considered in this work, it is necessary to modify the algorithm in order to work with noisy data.

Let us consider the relaxed version of (4.21) where we introduce a new variable ξ ∈ C m to be able able to capture the case when range(A) ∩ M b = ∅.

min (y,z,ξ)∈C m ×C m ×C m 1 2 ξ 2 s.t. y -z = ξ y ∈ range(A), z ∈ M b (4.22)
Let us define the update formulas for the Altenating Directions Method of Multipliers which is applied to (4.22).

Proposition 10 Let us consider the minimization problem (4.22). Then the Altenating Directions Method of Multipliers update formulas are

z (k+1) = b exp i arg(y (k) -ξ (k) + λ (k) ) , y (k+1) = AA † (z (k+1) + ξ (k) -λ (k) ), ξ (k+1) = ρ 1 + ρ y (k+1) -z (k+1) + λ (k) , λ (k+1) = λ (k) + y (k+1) -z (k+1) -ξ (k+1) ,
where

A † = (A * A) -1 A * is a pseudo-inverse of matrix A ∈ C m×n .
Proof. The scaled augmented Lagrangian writes as

L ρ (y, z, ξ, λ) = 1 2 ξ 2 + Re{ λ * (y -z -ξ)} + ρ 2 y -z -ξ 2 = 1 2 ξ 2 + ρ 2 y -z -ξ + λ 2 - ρ 2 λ 2 ,
where λ = λ ρ is a scaled dual variable. Then, we conclude that 1. The update formula for z with fixed y = y (k) , ξ = ξ (k) , and λ = λ (k) writes as

z (k+1) = argmin z∈M b ρ 2 y (k) -z -ξ (k) + λ (k) 2 = b exp i • argmin θ∈R m ρ 2 y (k) -b e iθ -ξ (k) + λ (k) = b exp i • arg(y (k) -ξ (k) + λ (k) ) .
2. The update formula for y with fixed z = z (k+1) , ξ = ξ (k) , and λ = λ (k) writes as

y (k+1) = argmin y∈range(A) ρ 2 y -z (k+1) -ξ (k) + λ (k) 2 = A argmin x∈C n ρ 2 Ax -z (k+1) -ξ (k) + λ (k) 2 = AA † (z (k+1) + ξ (k) -λ (k) ).
3. The update formula for ξ with fixed z = z (k+1) , y = y (k+1) and λ = λ (k) is obtained from

∂ ξ L ρ (y, z, ξ, λ) = ξ -ρ(y -z -ξ + λ) = 0,
from where we have

ξ (k+1) = ρ 1 + ρ y (k+1) -z (k+1) + λ (k) .
4. The update formula for λ with fixed z = z (k+1) , y = y (k+1) , and ξ = ξ (k+1) writes as

λ (k+1) = λ (k) + y (k+1) -z (k+1) -ξ (k+1) .
Observation 5 Let us consider the update formulas in Proposition 10. If λ (0) ∈ ker(A * ) then

λ (k) , ξ (k) ∈ ker(A * ) for all k ∈ N\{0}.
Furthermore, ρλ (k) = ξ (k) , for all k ∈ N\{0}.

Proof. First, note that

λ (k+1) = λ (k) + y (k+1) -z (k+1) -ξ (k+1) = λ (k) + y (k+1) -z (k+1) - ρ 1 + ρ y (k+1) -z (k+1) + λ (k) = 1 1 + ρ y (k+1) -z (k+1) + λ (k) = 1 ρ ξ (k+1) .
Then,

λ (k+1) = λ (k) + AA † (z (k+1) + ρλ (k) -λ (k) ) -z (k+1) -ρλ (k+1) from where λ (k+1) = 1 1 + ρ λ (k) + AA † z (k+1) + (ρ -1)λ (k) -z (k+1) = 1 1 + ρ λ (k) + (ρ -1)AA † λ (k) + (AA † -I)z (k+1)
Now, let λ (k) ∈ ker(A * ), from previous equations it is clear that λ (k+1) ∈ ker(A * ) since AA † -I is a projection matrix on ker(A * ). And thus, from the first equation we obtain that ξ (k+1) ∈ ker(A * ) also. 2

The Observation 5 is used for simplification the update formula for y (k+1) and for elimination variable ξ.

The main question for Algorithm 9 is how to select ρ properly. Before going further, let us provide numerical simulations for ρ = 0.5. On Figure 4.9 it can be seen that the problem with convergence for noisy case is fixed, where the oscillations during the convergence for the problems without noise in the model still exists but with lower frequency. It is also important to note that not all problems are solved by Algorithm 9 when ρ = 0.5. It appears that ρ = 0 is useful when an approximation x (k) is far from the solution, but then it must be changed to some positive value.

Let us provide the following observations which helps to understand what strategy for ρ is efficient.

Observation 6 If ρ = 1, in Algorithm 9, then it reduces to Algorithm 7 starting from the second iteration.

Proof. Follows from the second property of Observation 5.

2

Algorithm 9: Relaxed Alternating Directions Method of Multipliers. Input: Initial point x (0) ∈ C n and parameters ρ > 0, λ (0) , ξ (0) ∈ C m . Output: Approximate solution x ∈ C n of (1.2).

1 k = 0, y 0 = Ax (0) ; 2 Loop 3 z (k+1) = b exp i arg(y (k) + (1 -ρ)λ (k) ) ; 4 x (k+1) = A † z (k+1) ; 5 y (k+1) = Ax (k+1) ; 6 if stopping test is satisfied then 7 break 8 λ (k+1) = 1 1+ρ λ (k) + y (k+1) -z (k+1) ; 9 k = k + 1;
y -z ∈ ker(A * ) ∩ {αz, α ∈ (-ρ, ∞) m }.
Proof. Let (y, z, ξ, λ) ∈ C m × C m × C m × C m be a limit point of a sequence defined by update formulas in Proposition 10 exists and λ (0) ∈ ker(A * ), ρ > 0. Then, from Observation 5 we obtain that λ = 1 1+ρ λ + (AA † -I)z , y = AA † z, and ρλ = ξ which gives that y -z = ρλ.

Also, this can be equivalently written as AA † z -z = ξ. The last equation is used to obtain

z = b exp i • arg(AA † z -ξ + λ) = b exp i • arg(AA † z -AA † z + z + λ) = b exp i • arg(z + λ) , which gives that λ ∈ {αz, α ∈ (-1, ∞) m } and thus y -z ∈ range(A) ⊥ ∩ {αz, α ∈ (-ρ, ∞) m }. 2
Observation 8 Let us consider a recurrent form x (k+1) = 1 2 (x k +∆) for some starting x (0) ∈ R n and fixed ∆ ∈ R n . Then lim k→∞ x (k) = ∆.

Proof. The proof follows from the explicit form

x (k+1) = 1 2 k+1 x (0) + ∆ + ∆ k i=1 1 2 i . 2
Proposition 12 The Algorithm 9 always converges to a limit point (y, z, ξ, λ)

∈ C m × C m × C m × C m if ρ = 1.
Proof. The result for (y, z) ∈ C m × C m follows from Observation 6 and Proposition 7. The result for (λ, ξ) ∈ C m × C m follows from Observation 8 where ∆ = y -z and x (k+1) = λ (k+1) . 2

The formulation (4.22) is useful when range(A) ∩ M b = ∅, which indicates the presence of noise in data. The Observation 6 and Observation 7 give us an idea about the link between a regularization parameter ρ and the ability of algorithm to converge. However, it is not clear how to choose ρ properly at each k-th iteration of Algorithm 9. In the following part, we try to answer this question.

Proposition 13 Let ξ = AA † z -z be such that ξ ∈ {α z, α ∈ (-ρ, +∞) m } for z ∈ C m and ρ > 0. Then ξ -α z = 0 for

α = |z| -2 Re{AA † z z} -1. (4.23)
Proof. The optimal α can be found from the following optimization problem

min α∈R m f (α) := 1 2 ξ -α z 2
The constrains α ∈ (-ρ, +∞) m are not taken into account since it is given that a solution exists and it is unique due to strict convexity of f . Then, from the first order optimality conditions we obtain

∇f (α) = ∇ 1 2 ξ 2 -Re{ξ * (α z)} + 1 2 α z 2 = -Re{ ξ z} + α |z| 2 = 0.
Replacing ξ by AA † z -z the proof follows. 2

From Proposition 13 we can conclude that for any z ∈ C m the optimal α writes as It can be seen that ρ changes the solution set for z. For example, if ρ is close to zero, then vector z and its projection on range(A) must be almost colinear to satisfy AA † z -z -α z = 0.

α = max{-ρ 1 m , |z| -2 Re{AA † z z} -1}. ( 4
To have a clear picture in mind, let us visualize the set of possible projections for different ρ in the case of dimension 1. It can be seen on Figure 4.10 that when ρ = 0 the set of possible projections reduces to one point where the case for ρ = 1 indicates all possible projections in the set. The reason why we have such form of the set of possible projections can be seen on Figure 4.11 where the projection point is denoted by red color for many different linear spaces, which are visualized by gray dashed lines. This visualization brings an inside about how ρ affects the set of possible projections. However, at the same time, it gives an idea of how ρ can be properly selected during the iterations of the Algorithm 9. It is clear that at the begining of the iteration process it is reasonable to set ρ = 0 to guide iterations to the "perfect" point even if it does not exists. Closer algorithm converge to it -bigger ρ must be to be sure that the limit point exists. This kind of information can be obtained from (4.24). It is clear that to remove constrains defined by ρ, we need to have that -ρ = min{α 1 , . . . , α m , } for α i = |z i | -2 Re{(AA † z) i z i } -1 for i ∈ {1, . . . , m}. Equivalently, we have that ρ = max{-α 1 , . . . , -α m , }. The value of ρ tells us how far we are from solution. Thus, following the logic explained above, we can use it to construct the formula for ρ in the following form ρ = 1 -min{1, max{-α 1 , . . . , -α m }}.

(4.25)

The numerical study in Section 4.8 reveals the efficiency of this choice for ρ with a sufficient level of noise. However, when the noise dominates, the behavior near solution starts being unstable. For this goal, the empirical stabilization is proposed in the following form

If dist norm (y, z) < γ then ρ = 1, (4.26)
where γ ∈ [0, 1] depends on the level of noise in the data. Thus, the final form of ADMM is presented in Algorithm 10 where all useful propositions and observations are applied.

Algorithm 10: Adapted Relaxed Alternating Directions Method of Multipliers.

Input: Initial point x (0) ∈ C n , γ ∈ [0, 1]. Output: Approximate solution x ∈ C n of (1.2). 1 k = 0, ρ 0 = 0, ξ (0) = λ (0) = 0, y 0 = Ax (0) ; 2 Loop 3 z (k+1) = b exp(i arg(y (k) + (1 -ρ (k) )λ (k) )); 4 x (k+1) = A † z (k+1) ; 5 y (k+1) = Ax (k+1) ; 6 if stopping test is satisfied then 7 break 8 ρ (k+1) = 1 -min 1, max{|z (k+1) | -2 Re{y (k+1) z (k+1) } -1} ; 9 if dist norm (y (k+1) , z (k+1) ) < γ then 10 ρ (k+1) = 1 11 λ (k+1) = 1 1+ρ (k+1) λ (k) + y (k+1) -z (k+1) ; 12 k = k + 1; 13 return: x = x (k) .
Let us provide numerical simulations to verify the robustness and capability to solve all problems for Algorithm 10. On Figure 4.12 it can be observed that all problems almost without any oscillations can be solved by Algorithm 10 for the case when there is no noise in the model. However, there are still oscillations for the noisy case, when the algorithm reaches a limit of improvement of the approximation. That is why it is proposed to use γ parameter to switch to alternating projections when we are close to the solution. It was observed empirically that if γ = 2σ, then we can avoid the oscillations (Figure 4

Initialization methods

For nonconvex problems it is highly important to select a proper initial guess to avoid local minimums during iterations. There are several methods [START_REF] Gerchberg | Practical algorithm for determination of phase from image and diffraction plane pictures[END_REF][START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF][START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF][START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF][START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF][START_REF] Chen | Phase retrieval with one or two diffraction patterns by alternating projections of the null vector[END_REF] which all are from the family of spectral algorithms because they are based on the computation of the largest eigenvalue of a matrix

Y = 1 m m i=1 f (b i )a i a * i
, where b i ∈ R + is the i-th measurement and a i ∈ C n is the i-th row of a transmission matrix A ∈ C m×n such that b = |Ax|, and function f : R → R differs for each method. The proper choice of f can reduce the required m to estimate x [START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF][START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF][START_REF] Chen | Phase retrieval with one or two diffraction patterns by alternating projections of the null vector[END_REF]. Since the goal of this work is to create a fast and robust algorithm to solve a phase correction problem, in this section we consider just two of initialization strategies to use in further numerical experiments instead of taking into account all of them. Following the published results in [START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF][START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF], we decide to select Wirtinger flow (WF) initialization method [START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF] as a baseline algorithm, and Gao & Xu (GX) initialization algorithm [START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF] which is proved in the original paper to have better properties than WF. Both WF and GX initialization algorithms requires computation of the leading eigenvector of matrix Y . The power method can efficiently solve this problem with a selected accuracy. Despite the fact that it can be easily found in literature, we put it in Algorithm 11 to make this document self-contained.

Algorithm 11: Power method.

Input: Matrix M ∈ C n×n , tolerance ε > 0.

Output: The largest eigenvalue λ ∈ C and eigenvector v ∈ C n .

1 Set v = 1 n n , λ = v * M v v 2 . 2 Loop 3 Set v = M v M v . 4 Set λ = v * M v v 2 . 5 if |λ -λ| < ε then 6 break 7 Set λ = λ . 8 Set λ = λ ; 9 return λ, v.
The goal of this section is not only to present the algorithms but also give the idea behind them. The principle of spectral methods is to use a transmission matrix A ∈ C m×n and a vector of measurements b ∈ R m + such that b = |Ax| to approximate x ∈ C n . The Algorithm 12 is based on the fact if A is Gaussian then by the law of large numbers we have

1 m m i=1 b 2 i a i a * i ≈ x 2 I + xx * . (4.27)
It can be seen that any leading eigenvector of x 2 I + xx * is of the form λx for λ ∈ R. Thus for sufficiently large m, solution x can be recovered up to a sign and a global phase shift. The normalization coefficient λ approximates x but has no effect for a class of projection algorithms that is considered in this work. The reason comes from the update formula in Algorithm 7 that uses arg(Ax) expression, which is clearly invariant to multiplication by a real positive constant.

Algorithm 12: Wirtinger Flow: Initialization [7, Algorithm 1].

Input: Vector of measurements b ∈ R m + , transmission matrix A ∈ C m×n , Power method tolerance ε > 0.

Output: Initial guess x (0) ∈ C n . 1 Set Y = 1 m m i=1 b 2 i a i• a * i•
, where a i• is the i-th row of A.

2 Find the larges eigenvalue v ∈ C n of Y using Algorithm 11 with tolerance ε.

3 Set λ = n m i=1 b 2 i m i=1 a i• 2 . 4 Set x (0) = λv.
The main result about WF initialization is presented in Theorem 8 which claims that O(n log n) measurements are required to obtain an efficient initialization with high probability.

-1 m i=1 b 2 i , 1 obeys dist(x (0) , x) ≤ 1 8 x ,
with probability at least 1 -10e -γn -8/n 2 (γ is a fixed positive numerical constant).

The idea of Gao & Xu initialization algorithm is similar to WF but with f (x) =1 2 -exp(x/λ 2 ) which improves the quality of initial guess estimation. In the same way that we used in Lemma 1 let us show the motivation of such f .

Lemma 2 [13, Lemma 5.2]. Let a = Re[a] + i Im[a] such that Re[a], Im[a] ∼ N 0, 1 2 , then for any fixed vector x ∈ C n it holds that E 1 2 -exp |a * x| 2 /λ 2 aa * = xx * 4 x 2 .
The main result about GX initialization is presented in Theorem 9 which claims that O(n) measurements are required to obtain an efficient initialization with high probability. Input: Vector of measurements b ∈ R m + , transmission matrix A ∈ C m×n , Power method tolerance ε > 0.

Output: Initial guess x (0) ∈ C n . 1 Set λ 2 = m i=1 b i m . 2 Set Y = 1 m m i=1 1 2 -exp b i /λ 2 a i• a * i•
, where a i• is the i-th row of A.

3 Find the larges eigenvalue v ∈ C n of Y using Algorithm 11 with tolerance ε.

4 Set x (0) = λv.

Theorem 9 [13, Theorem 2.1] Let x ∈ C n be an arbitrary vector and b = |Ax| ∈ R m + be m samples. Then for any θ > 0 there exists constant C θ such that for m ≥ C θ n the computed by Algorithm 13 estimate x (0) normalized to have squared Euclidean norm equal to m

-1 m i=1 b 2 i , obeys dist(x (0) , x) ≤ √ 3θ x ,
with probability at least 1 -4e -c θ n (c θ is a fixed positive numerical constant).

Let us now provide a few time profiles for the phase retrieval methods that we discussed above which use different initialization strategies: random (RI), Wirtinger (WI) and Gao & Xu (GXI) initializers. It appears that they have different impact on the convergence of the considered phase retrieval algorithms. For example, on Figure 4.14a it can be seen that both strategies WI and GXI increases the number of solved problems in comparison with a random starting point. However, both of them have approximately the same effect on the speed of convergence for the gradient descent algorithm with a secant line search. The secant line search was selected since it is faster than backtracking line search, which can be seen on Figure 4.3. This picture is not the same for other methods. For instance, on Figure 4.14b it can be seen that GXI initialization allows to solve more problems than both WI and RI, where at the same time for half of the problems WI allows to converge faster. Alternating projections algorithm reveals that GXI initialization improves both the speed of convergence and the number of solved problems in comparison with Wirtinger and random initialization strategies. The most interesting behavior can be observed on Figure 4.14d. It reveals that the strategy which allows both to obtain the fastest convergence and to solve all problems is a random initialization strategy. Recall, that WI and GXI requires additional computational time to find starting point. More precisely, it needs to compute a matrix Y and to find its leading eigenvector by means of power method. It reveals that this effort is not relevant for ADMM (Algorithm 10) and it is more efficient from the computational speed point of view to generate randomly an initial point.

The experiments were done for n = 32, m = 128, without noise σ = 0, where 100 problems were generated for this aim. The metric dist norm is used to measure distance between a solution and an approximation. Also, the stopping tolerance ε in Algorithm 11 was set to 0.001 in order to approximate the leading eigenvector.

Numerical results

In this section, we present numerical results for different optimization algorithms to solve a set of phase retrieval problems in order to find the fastest and the most robust algorithm.

Following the previous research, we exclude Gauss-Newton method from a list of comparison. The reason can be clearly seen on Figure 4.14, where Gauss-Newton method can be up to 10 times slower than other algorithms. Also, we use different initialization strategies for different algorithms. On Figure 4.14 it can be clearly seen that the best initialization approach for gradient descent is Wirtinger flow, alternating projections algorithm performs faster with Gao & Xu initialization, and ADMM does not require any initialization algorithms and uses a random complex vector as a starting point.

We present time and iterations profiles for these three algorithms with corespondent initialization strategies, for different number of beams n ∈ {8, 32, 128} and measurements m = 4n, for two levels of noise σ = 0 and σ = 0.1. The stopping tolerance of metric dist norm for σ = 0 is 0.001, for σ = 0.1 is 0. These profiles reveals that ADMM (Algorithm 10) in all of the cases performs faster than alternating projections (Algorithm 7) and gradient descent (Algorithm 5) algorithms. In addition, with increasing n this difference in speed increases. Now, let us provide heat maps (explained in Section 2.6) for these algorithms to observe the solving capabilities for different n and m, with and without noise. It can be clearly seen (Figure 4.21 -Figure 4.26) that alternating projections and gradient descent algorithms requires m ≥ 4n measurements to solve a phase retrieval problem for both cases: σ = 0 and σ = 0.1. At the same time, ADMM algorithm can solve more that 50% of all problems even when m = 3n, which can be observed on Thus, we can conclude that ADMM (Algorithm 10) appears to be faster and more robust (with a properly selected value of γ) than other algorithms considered in this work, which potentially gives a real improvement of the phase correction algorithm considered in the next chapter.

Chapter 5 Phase corrections by means of optimization methods

In this chapter, we consider a phase correction problem to solve using a set of optimization algorithms, which aim to find a solution of (1.2). In Section 1.3, we briefly described the problem and now, let us give a formal definition.

Let us consider a vector x (k) ∈ C n representing n input lasers whose vector of phases arg(x (k) ) is unknown, but with a known vector of amplitudes |x (k) |. The output is given by a vector of m intensity measurements b (k) ∈ R m + for m > n. The link between x (k) and b (k) can be modeled by means of a transmission matrix A ∈ C m×n such that b (k) = |Ax (k) | where A represents a scattering device and operation | • | simulates the measurement process. Let x ∈ C n be a given target signal. The goal is to apply a phase correction to the vector x (k) to obtain a new vector x (k+1) such that |x

(k+1) | = |x (k) | and |Ax (k+1) | = |A x|. Suppose that x (k) is a solution of the equation |Ax| = b (k) , where b (k) = |Ax (k) |, (5.1)
such that arg( x (k) ) = arg(x (k) ) up to a constant. The current phase values are then modified by means of the phase modulation to produce a new vector x (k+1) such that arg(x (k+1) ) = arg(x (k) ) -arg( x (k) ) + arg( x).

Therefore, for some constant c ∈ R, we will have

arg(x (k+1) ) = arg( x) + c.
The solution of (5.1) can be done by a phase retrieval algorithm (for instance, the alternating projection). But the computation of a nearly exact solution of this equation is not relevant in our context. The first reason is that the computational cost can be too large. The second reason is that the convergence to a global solution is not always guaranteed. The third reason comes from robustness. Due to noisy data, it is useless to compute a nearly exact solution. It is better to perform some few iterations of a phase retrieval algorithm to get an approximate solution of (5.1), to modify the input sources by phase modulation, then to perform new magnitude measurements and then to loop. This process is formulated in Algorithm 14. In practice, this loop is continuously applied and never stop because of phase deviations.

It is worth noting that the steps 2 and 4 are not numerical implementations and that the value of arg(x (k) ) is never explicitly known. In Step 2, the vector b (k) follows from physical measurements by means of a scattering device with photodetectors.

Step 4 corresponds to phase modulations of the input sources. For a detailed experimental setting and true experiments, see [START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF]. However, in our numerical experiments the phases are given by the arguments of a true vector x (k) and the magnitude measurements are simulated by calculating |Ax (k) |.

Numerical results

In this section, we present numerical results for different phase retrieval algorithms applied to the problem of a phase correction. At first, let us select the best initialization strategy for starting point computation in the same way, as we did on Figure 4.14. The reason why the results on Figure 4.14 are not used directly in this section is that now we solve a phase correction problem, but not a phase retrieval, and the metric to measure distance between a solution and an approximation is q norm instead of dist norm . Thus, from Figure 5.1 we conclude that:

1. Gauss-Newton is extremely slower than other methods (Figure 5.1b) and will not be considered for further comparison. The experiments were done for n = 32, m = 128, without noise σ = 0, where 100 problems were generated for this aim. The metric q norm is used to measure distance between a solution and an approximation. Then, from Figure 5.2 -Figure 5.7 we conclude that ADMM algorithm, which is used in step 3 of Algorithm 14 has the highest speed in solving a phase correction problem in both cases: with noise σ = 0.1 and without it (σ = 0).

Phase corrections by means of neural networks

In this chapter, the Algorithm 14 is considered where a neural network model is used in step 3 to compute an approximation x (k) . Two kind of results were obtained in this domain. The first is the way to train a neural network by the quasi-reinforcement learning algorithm (Algorithm 2) and use it in a phase correction loop where a target x ∈ C n is fixed. It is shown that this model is highly scalable and needs extremely low time to compute x (k) ∈ C n . However, the practical application requires the possibility to dynamically change x which is impossible in this context. To overcome this difficulty, the second result is presented where it is revealed how to build and train a target adaptive neural network. The advantages and disadvantages are discussed also.

Research path

The form of the networks that are considered in this chapter have a simple linear nature, which may be confusing at the first glance since the problem that they solve is complex and nonlinear. However, their simplicity arises from an among of numerical and experimental tests of different network architectures, and several learning approaches.

At the beginning, the idea was to train a deep network with many layers and different activation functions to solve directly a phase retrieval problem, which would have led to a phase correction in one step. However, this idea works just for a small number of beams (up to 6). That is why it was proposed to integrate a network in a phase correction loop in order to train to perform a set of corrections instead of just one. Because of a high similarity with the reinforcement learning framework we used its terminology to explain our idea.

Due to the limitation in a phase correction speed, the network must be lightweight and simultaneously capture nonlinear dependencies. The first idea was to use a neural network with two layers and ReLU activation function between them to compute phase correction values. However, the quality of this model was not appropriate for our goals.

In addition, it was noted that instead of predicting phase values, it is better to predict real and imaginary parts of a complex number and then apply arg(•) function to retrieve phases. This has a logical reason because instead of predicting fixed values ϕ ∈ [0, 2π], we aim that the relation of imaginary and real parts be fixed. We can say that this approach has more degree of freedom in some sense.

A prediction of real and imaginary parts and a usage of arg(•) function at the output gives not only an additional degree freedom, it also naturally introduce a nonlinearity to the network. Note, that even if the model is linear, the usage of arg(•) can be considered as a nonlinear activation function. It was observed numerically, that there is no need to use a multilayer neural network with classical activation functions because of the presence of arg(•). That is why the simple network form (6.1) is proposed. S t is the vector of intensities b (k) and the reward R can be computed by q norm between x (k) and x (k) . The adaptation of the Algorithm 2 is presented in Algorithm 15, where each step in Algorithm 15 is connected to some step in Algorithm 2. Let us explain each step of Algorithm 15 in details. To simplify the notations for explanations of steps 2-7 let us consider that k is fixed and remove it from notations, and consider that j ∈ {1, . . . , N } where it is not explicitly specified.

In step 1 the initial parameters N , K max , I max , x, and W are initialized. Let us give the information about all of them. The size of batch N affects the speed of convergence of the learning algorithm. A small value like 16, 32 or 64 can produce low informative vector of the ascent direction g in step 5. A large value like 10 4 or 10 5 gives to much information for g. The problem is that the function to optimize (in our case it is f (W ) := 1 N N j=1 q norm e iΦ j , X j ) is not convex which means that if g is computed with big N then with high probability we converge to one of local minimums. That is why there is an empirical trade off between little and highly informative g. It was observed that for N = 4096 this trade off is satisfied (see Figure 6.1). The second parameter to explain is a maximum number of simulated corrections K max . The idea of the Algorithm 15 is to find parameters of (6.1) such that

1 N N j=1 q norm e iΦ (Kmax) j , X (Kmax) j ≈ 0. ( 6.3)
Let K max = 1, then the Algorithm 15 learns the model to find a nearly exact solution of a phase retrieval problem (5.2). However, it was observed empirically that for large n the problem (5.2) is difficult to solve (see Figure 6.

2).

That is why it is proposed to update the parameters of model (6.1) at each of the k-th correction for K max > 1 in order to improve the approximation of a solution for (5.2). The exact value of K max depends on both number of beams n and number of measurements m. The dependence is such that if n increases, required K max increase also. However, if m is increased for constant n, then K max decreases. It was observed numerically that K max = 9 is sufficient for training up to n = 128 and m = 4n. (see Figure 6.3). The number of learning iterations I max is also an empirical value which is set to I max = 1000 in the experiments. However, if n is large and it is observed that the value of q norm still decreasing, then I max must be increased. The Algorithm 15 also requires a target x as an input parameter. This point is crucial for understanding because just considering a fixed target the learning process converges. The trained model (6.1) can be seen as an approximation of the inverse of x → |Ax| in the neighborhood of x which is formulated in the Conjecture 1. This leads to the fact that if the target is changed, then the model (6.1) must be retrained. The initial parameters W are distributed uniformly on the range [-l, l] where l = 6/(m + n). This is a standard way to initialize the parameters of a model which is called Glorot uniform initializer or Xavier uniform initializer [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF].

In step 2 the batch of random signals X = {x j } N j=1 ⊂ C N ×n is generated where |x j | = | x| and arg x jk ∼ Uniform(-π, π).

In step 3 the vectors of intensity measurements b j ∈ R m + for each x j are computed. This step can be performed either numerically or experimentally. For numerical measurements, a model like (1.2) is required. The advantage of this approach is a computational speed. The disadvantage is that the model must be built. For experimental measurements, a fast phase modulation physical device is required. The advantage is that we do not need a model. The disadvantage is that in practice to generate batch of signals X ∈ C N ×n a Spatial Light Modulator (SLM) is used. This devise has low phase modulation speed around 500 milliseconds, which significantly slows down the learning process.

In step 4 the vectors of phase corrections ϕ j ∈ [-π, π] n are computed for each vector of measurements b j ∈ R m + . There is a link between a vector of phases ϕ obtained in step 4 of the Algorithm 15 and x in step 3 of the Algorithm 14 is such that ϕ = arg( x).

In step 5 the gradients are computed separately with respect to real and imaginary parts of parameters W ∈ C m×n . Despite the fact that the matrix is complex and potential application of Wirtinger calculus (see Section 2.3) to compute a complex-valued gradient is possible, the real-valued gradient is used instead. The motivation is that it was observed empirically that the convergence of the learning process is better with a real-valued gradient (see Figure 6.5)

The explicit formula of gradient is not given because it is computed by means of automatic differentiation algorithm implemented in Tensorflow library.

In step 6 the parameters W are updated using a real-valued gradient and following the update rule defined in [18, Algorithm 1] with default parameters α = 0.001, β 1 = 0.9, β 2 = 0.999. In step 7 the phase corrections ϕ j are applied for each x j . In (6.2) the subtraction is considered between arg( x) ∈ [-π, π] 1×n and rows of Φ arg

( x) -Φ =   arg( x) -ϕ 1• . . . arg( x) -ϕ N •   .
In step 8 it is checked if the maximum number of corrections K max is reached. If it is true, then a new batch X ∈ C N ×n is generated.

Also, there are several important properties of (6.1) which is necessary to mention. Let us show that the function (6.1) is invariant to changes of the absolute values of the measured signal. For some applications, it is possible that the transmission matrix must be normalized as it is defined in Definition 5 to remove the bias which comes from the physical system. Thus, it is important to be normalization invariant for a phase correction model. 2

By the Proposition 15 we conclude that the model (6.1) is invariant to a transmission matrix normalization.

To obtain a full picture regarding the capabilities of TFNN, several additional information slices are presented in Figure 6.6.

It was numerically observed that in order to achieve a sufficiently low q norm , say q norm < 0.04, there is a minimal required ratio m/n for different n, which equals m = 4n. When the beam count varies from 4 to 128, the required m/n ratio increases from 2 to 12. Thus, it is important to show a minimal required ratio m/n for different n to achieve a sufficiently low q norm . Different TFNNs were trained for the various number of beams n ∈ {4, 8, 16, 32, 64, 128} and the different ratios between the number of measurements and the number of beams m/n ∈ {2, 4, 6, 8, 10, 12}. The minimal achievable q norm was recorded and visualized as a heat map in Figure 6.6a, with the corresponding relative training time shown in Figure 6.9b. The minimal achievable q norm is obtained by solving 1000 phase correction problems with fixed targets for each combination of n and m/n, and computing average of the q norm at the last correction. Also, an average q norm during corrections is visualized in Figure 6.7, which gives the information about the number of corrections required to achieve a target.

Target adaptive neural network

In the previous section it is shown how to train a neural network to use in the phase correction loop with a fixed target. However, it is impossible to retrain the model (6.1) fast when this the number of targets in a generated batch. The parameter M represents the number of initial signals generated for each of N targets. The exact values that are used in computations are N = 1024, M = 256. Matrix U ∈ C mn×n is initialized such that the entries of real and imaginary parts are distributed by a standard normal law.

In step 2 the sets of target signals X ∈ C N ×n and initial signals X ∈ C N ×M ×n are generated such that | X| = |X j | = c j for some constant c j > 0, arg(X jpq ), arg( X jq ) are distributed uniformly on range [-π, π] for q ∈ {1, • • • , n}.

In step 3 matrices W j ∈ C n×m are computed for each target signal x j ∈ C n by means of model (6.4).

In step 4 the measurements B j ∈ R M ×m for each X j ∈ C M ×n are performed. The same as in the explanation of Algorithm 15, this step can be performed either numerically or experimentally.

In step 5 the phase corrections Φ j ∈ [-π, π] M ×n are computed for each B j ∈ R m + by means of matrices U j ∈ C n×m .

In step 6 the gradients are computed. Note, that Wirtinger calculus is used to get g in contrast to the Algorithm 15. It was observed numerically that the real-valued gradient with respect to real and imaginary parts of U is less informative for optimization (see Figure 6 In step 7 the matrix W is updated. Note, that the default parameters of Adam update rule [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] are changed to α = β 1 = β 2 = 0.1. It was observed empirically that the optimization with complex-valued gradient is faster with these parameters.

In step 8 phases of each X j ∈ C M ×n are corrected with Φ j ∈ [-π, π] M ×n . To obtain a full picture regarding the capabilities of TANN, several additional information slices are presented in Figure 6.9. It was numerically observed that in order to achieve a sufficiently low q norm , say q norm < 0.04, there is a minimal required ratio m/n for different n. When the beam count varies from 4 to 20, the required m/n ratio increases from 4 to 12. Thus, it is important to show a minimal required ratio m/n for different n to achieve a sufficiently low q norm . Different TANNs were trained for the various number of beams n ∈ {4, 6, 8, 10, 12, 14, 16, 20} and the different ratios between the number of measurements and the number of beams m/n ∈ {2, 4, 6, 8, 10, 12, 14, 16, 18, 20}. The minimal achievable average q norm Chapter 7

Experimental results

In this chapter we compare the efficiency and robustness of the phase retrieval methods used in the experimental phase correction loop. In this work, we developed two kinds of such methods: Adaptive Relaxed Alternating Directions Method of Multipliers (Algorithm 10) to solve a phase retrieval problem and Target Adaptive Neural Network (6.4) (learned by Algorithm 16) to solve directly a phase correction problem. Both of them were applied to correct n = 16 laser beams with 4n measurements for ADMM and 12n measurements for TANN. Both of them were compared with a baseline alternating projections algorithm (Algorithm 7). The initialization strategy was not used when ADMM and alternating projections algorithms were compared, which means that the same random starting point was used instead. The reason is that it was observed on Figure 5.1 that both algorithm have the same kind of profit from the spectral initialization.

To compare the algorithms the profiling method (Section 2.5) is used with modifications. In the initial approach, we generate a set of problems (the set of transmission matrices) that are solved by selected methods. However, in practice, it is too complicated to perform such computations since around 10 minutes is required to collect experimental data X ∈ C N ×n and B ∈ R N ×m + for N = 400 to compute a transmission matrix A. That is why we generate different initial state of lasers instead and then apply the phase correction loop for a fixed transmission matrix A ∈ C m×n .

The transmission matrix was computed by the Algorithm 17 using unbiased data B -β, where β was computed by the Algorithm 18. The same matrix was used to train TANN in Algorithm 16.

The experimental profiles are presented on Figure 7.1 and Figure 7.3 reveal that both algorithms performs faster that the baseline algorithm and both can solve all problems. ADMM was applied for m = 4n measurements with stopping tolerance 0.01 of q norm metric and 15 iterations maximum and solve 80% of all problems in 0.5 milliseconds where, at the same time, the alternating projections with stopping tolerance 0.001 of q norm metric solves 15%. The stopping tolerances were selected empirically in order to solve a phase correction problem as fast as possible. It was observed that these values are different for different algorithms. The speed of TANN is incomparably greater than for the alternating projections algorithm, which is obvious since one correction requires a single matrix multiplication. As it was shown on Figure 6.9a, TANN requires more measurements than optimization algorithms, which is a disadvantage. However, when this amount is available, than the highest speed of corrections can be achieved. The experimental setup that was built by PhD student Alexandre Boju and postdoctoral researcher Geoffrey Maulion is presented on Figure 7.5. The details about the experimental setup are described in the thesis of Alexandre Boju.

Conclusion and perspectives

Our main contributions in this work were first to develop the ADMM algorithm (Chapter 5) and to show that this algorithm can outperform the alternating projections method initially used to perform phase corrections in the context of physical experiments in the photonics laboratory (Chapter 7). In particular, we have observed and shown that the choice of the penalty parameter is sensitive for the efficiency of the algorithm and we have proposed adaptive updating rules that we believe are relevant. Our second main contribution was the development "from scratch" of a neural network approach (Chapter 6) for the phases control. The first algorithm (ADMM) has been developed to solve the phase correction problem, but it can also be used to solve the phase retrieval problem (Algorithm 10). The updated formulas of an algorithm and its speed of convergence depends on the reformulation of the original system of nonlinear equations |Ax| = b. The proper reformulation was found, which allowed us to achieve a higher speed in solving a phase correction task. It was also shown how to select a regularization parameter ρ of the augmented Lagrangian to speed up the algorithm and, at the same time, maintain the ability to converge to a global solution frequently. To summarize, there are the following advantages and disadvantages of ADMM in comparison with the baseline Alternating Projections (AP) algorithm. The comparison is valid for a number of beams n up to 128 and for fixed m = 4n. Advantages:

1. ADMM requires less time and iterations to solve both a phase retrieval and a phase correction problems.

2. ADMM converges to a global solution of a phase retrieval problem with higher probability than AP.

Disadvantages:

1. ADMM requires a proper selection of the parameter γ which is used to switch to the AP algorithm when it is near a solution of a phase retrieval problem. This must be used when there is noise in a transmission matrix.

The second algorithm solves a phase correction problem only and is based on the neural network. The main contribution is not the form of a model but the way of its training : Algorithm 15 for a fixed target phase and Algorithm 16 for a random target phase. Since a target phase changes with high frequency in the physical application we will consider model (6.4) only in this conclusion. Let us reveal the advantages and disadvantages of this approach. Advantages:

1. Because of a simple form of the neural network (6.4), the corrections are performed with high speed. If we compare the profiles for experimental data (Figure 7.3), it will be clear that TANN is approximately 15 times faster than AP, which is the biggest advantage of this approach.

2. Figure 7.3 reveals that the neural network model is robust to noise.

Disadvantages:

1. To train TANN model (6.4) we require much more measurements than for the optimization algorithms. It can be clearly seen on Figure 6.9a.

2. TANN model has memory limitations. For example, to train a model for 20 beams, we require GPU with 6GB. Taking into account that raises is quadratic, this algorithm can not be applied for a huge number of beams.

To conclude, both methods perform faster than the alternating projection algorithm. However, there are limitations in memory for TANN, which makes this approach not applicable for n > 20, and additional parameter γ for the ADMM algorithm, which must be selected with respect to the noise level.

This conclusion leads us to the future perspectives of this work. Let us first specify them for the ADMM algorithm, and then for TANN.

From a practical point of view, ADMM requires the possibility to select the parameter γ. This can be easily done if the experimental data of signals X ∈ C N ×n and measurements B ∈ R m + is available. For this aim, the standard deviation σ of noise must be estimated statistically and then γ = 2σ (empirical result). However, the very last results from the physicists team give a possibility to compute a transmission matrix without X. This modification reduces the size of the laser system, which is important for real applications. That is why an approach for the computation of γ without X must be developed.

From a theoretical point of view, ADMM requires a convergence result. This must be done to fully understand the capabilities of this algorithm. However, since the original optimization problem is nonconvex and complex-valued, there are not a lot of available results in the literature, which makes this analysis very hard and time expensive. For instance, in [START_REF] Waldspurger | Phase retrieval with random gaussian sensing vectors by alternating projections[END_REF] there is a result for the convergence of the alternating projection algorithm. However, the assumptions that must be satisfied to achieve the convergence are not possible to meet in the real-life (millions of measurements), which is also mentioned by the author. That is why it is a really challenging problem to solve.

Concerning the neural network approach, the biggest problem is a memory limitation. The quadratic raise of the trainable parameters is a huge disadvantage. In this context, the complexvalued convolutional layer could be potentially useful. It is already implemented in the Tensorflow extension library, which is used in this work, however, the proper way of its application was not found. That is why it is reasonable to work in this direction in the future.

The results about the neural network approach were published in [START_REF] Shpakovych | On-demand phase control of a 7fiber amplifiers array with neural network and quasi-reinforcement learning[END_REF][START_REF] Shpakovych | Experimental phase control of a 100 laser beam array with quasi-reinforcement learning of a neural network in an error reduction loop[END_REF] and presented in the following conferences [START_REF] Shpakovych | Scalable and fast phase control of laser beam array with quasi-reinforcement learning of a neural network in an error reduction loop: simulation and experiment[END_REF][START_REF] Shpakovych | Combinaison cohérente de 100 faisceaux lasers grâce à un réseau de neurones artificiels appris par quasi-renforcement[END_REF][START_REF] Shpakovych | Phase control of a laser beam array, from optimization process to artificial intelligence[END_REF][START_REF] Shpakovych | Verrouillage et structuration dynamique des phases d'un réseau de 7 lasers par approche neuronale[END_REF][START_REF] Shpakovych | Optimization and machine learning algorithms applied to the phase control of an array of laser beams-renforcement[END_REF][START_REF] Shpakovych | Machine learning algorithm applied to the phase control of an array of laser beams[END_REF][START_REF] Shpakovych | Phase locking of fiber laser array using quasi-reinforcement learning, principle and experiments[END_REF], where [START_REF] Shpakovych | Combinaison cohérente de 100 faisceaux lasers grâce à un réseau de neurones artificiels appris par quasi-renforcement[END_REF][START_REF] Shpakovych | Phase control of a laser beam array, from optimization process to artificial intelligence[END_REF][START_REF] Shpakovych | Verrouillage et structuration dynamique des phases d'un réseau de 7 lasers par approche neuronale[END_REF][START_REF] Shpakovych | Phase locking of fiber laser array using quasi-reinforcement learning, principle and experiments[END_REF] were invited conferences.

The Python implementations of all algorithms that are presented in this work can be downloaded from the following links: inside the range [µ-z α s/ √ m, µ+z α s/ √ m] with probability P = 0.95, where z α = 1.96. Finally, assuming that all relative norms are drawn from the normal distribution with parameters μ and σ, we use a rule of 3-sigma which says that the sample ξ from this distribution is such that ξ ∈ [μ -3σ, μ + 3σ] with probability P = 0.9973. This computations assumes that most of the approximations a j are global solution and just some of them are local. This assumption is reasonable since it was shown numerically that with Gao & Xu initialization more than 90% of all problems can be solved by means the alternating projections algorithm. However, to be sure that the threshold is not too high we set a limitation on the level of 0.2 of relative norm. Then, the algorithm of finding A can be formulated in Algorithm 17. Find set J such that for all j ∈ J condition (A.2) is not satisfied; Note, that this algorithm gives no guarantee that rows a j of A for j ∈ {1, . . . , m} are global solutions of respective phase retrieval problems |Xa| = b •j . The are two reasons. First, we fix the maximal number of restarts R max , which is an empirical value and can be wrongly defined by user. Second, the maximal value of threshold 0.2 is also an empirical value, which means that we do not assume that the level of noise in measurements is greater than 20%. This can be not true in practice, however consequently it means that the measurements are too noisy and the model will have a low quality, which will lead to poor phase corrections. That is why, there must be warnings that the program can throw to inform user about problems with matrix computation and suggest possible solutions: either to increase R max or to inform about a high presence of noise.

When the transmission matrix is retrieved we can compute an error between experimental Then, to obtain an unbiased transmission matrix we first compute β by Algorithm 18, and then compute A by Algorithm 17 for X and B -β input data. Note, that to use this model, we require subtract β from experimental measurements each time when they are performed. The result of unbiasing can be observed on Figure A.1b. Also, the value of (A.3) can be compared. On In practice, it is possible that the transmission matrix must be normalized as it is defined in Definition 5 to remove the bias which comes from the physical system. The way to do that is presented below.

Figure 1 . 1 :

 11 Figure 1.1: Possible applications of a beam with a high brightness.

α

  j • e iϕ k → max {ϕ 1 ,...,ϕn}.

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: Phase correction to a zero target (MO: master oscillator, ϕ: phase modulator, A: fiber amplifier)

Figure 1 . 5 :

 15 Figure 1.5: Mathematical notations for physical values.

Figure 1 . 6 :

 16 Figure 1.6: Machine learning paradigm.

1. 3 .Figure 1 . 8 :

 318 Figure 1.8: One dimensional phase correction visualization.

Figure 2 .

 2 Figure 2.1: q(x, y)

Figure 2 . 2 :

 22 Figure 2.2: Profiles example.

Figure 2 . 3 :

 23 Figure 2.3: (a) Heat maps of the fraction of solved problems in grey scale and (b) its required average number of iterations.

Figure 2 . 4 :

 24 Figure 2.4: (a) Heat maps of the minimal achievable mean q norm in grey scale and (b) its required relative training time. The relative time on (b) is computed by dividing a learning time in seconds for each n and m/n by the minimal time to obtain GPU invariant information. The minimal time required by the GPU used in for this experiments was 1.3 s.

  ) where x ∈ F n is an input data and W ∈ F m×n , b ∈ F m are parameters to optimize. Sometimes, we will skip parameters W and b from the function arguments list and write fc(x) = W x + b to simplify notation. Note, that if m = 1, n = 2 and F = R, fully connected layer coincides with the function a(x, w) from Example 5 where w 0 = b.

Example 6

 6 Let ξ 1 ∼ N (-2, I 2 ) and ξ 2 ∼ N (2, I 2 ). Here I 2 is an identity matrix 2 × 2. Let us sample ξ 1 100 times and collect it into X 1 . Then we do the same for ξ 2 and X 2 . Thus X = [X 1 , X 2 ] ∈ R 200×2 is a set of 200 objects where each element x ∈ X was sampled either from N (-2, I 2 ) or from N (2, I 2 ) distributions. This process creates an artificial data from two classes. Then, let us prepare the responses. Let Y 1 ∈ 0 R 100 and Y 2 ∈ 1 R 100 which means that

Figure 3 . 1 :

 31 Figure 3.1: Generated objects X.

  is the set of responses. At this moment we created the training sample D = (X, Y ) of size 200. Let us define the parametric family of functions where we search the best one. Our model contains one fully connected layer fc : R 2

. 9 )

 9 Now, we split training sample on the train D train and test D test disjoint subsets and define X train , Y train are the objects and responses from D train , the same for X test and Y test . Then the problem of finding parameters (w, b) * can be formulated as the following optimization problem (w, b) * ∈ argmin (w,b)∈R 3 L(a(X train , w, b), Y train ).(3.10)

Figure 3 . 2 :

 32 Figure 3.2: Learning process.

Figure 3 . 3 :

 33 Figure 3.3: Separation hyperplane.

Figure 3 . 4 :

 34 Figure 3.4: The most common activation functions.

  a) Leaky ReLU function (α = 0.3). (b) ELU function (α = 1.0).

Figure 3 . 5 :

 35 Figure 3.5: ReLU function modifications.

Algorithm 1 :

 1 Supervised learning algorithm 1. Initialization: (a) Let w 0 ∈ W(F) be the starting value for trainable parameters. (b) Let a : X × W(F) → Y be a model to learn. (c) Let L : Y × Y → R is a loss function. (d) Let X train ⊂ X, Y train ⊂ Y are train objects and responses.

4 .

 4 Return loop: Set k = k + 1 and go to 2.

Figure 4 . 1 :

 41 Figure 4.1: (a) Metric dist norm between approximation x (k) computed by Algorithm 3 with line search Algorithm 4 and a solution for n = 8, m = 32, and σ = 0, and (b) for σ = 0.1.

9 if stopping test is satisfied then 10 break 11 k = k + 1 ;

 9111 Let us provide numerical simulations to observe how Algorithm 5 can solve phase retrieval problem (1.2) for n = 8 and m = 4n The results of numerical simulations are presented on Figure4.2, where the meaning of curves are the same as for Figure4.1. From Figure4

Figure 4 . 2 :

 42 Figure 4.2: (a) Metric dist norm between approximation x (k) computed by Algorithm 5 and a solution for n = 8, m = 32, and σ = 0, and (b) for σ = 0.1.

Figure 4

 4

Figure 4 . 3 :

 43 Figure 4.3: (a) Profile for solved problems by selected time for 100 problems of size n = 8, m = 32, and σ = 0, and (b) for σ = 0.1.

Algorithm 6 : 5 x 6 if stopping test is satisfied then 7 break 8

 65678 Gauss-Newton algorithm. Input: Initial point x 0 ∈ C n , regularization parameter δ > 0. Output: Approximate solution x ∈ C n of. 1 x = x 0 ; 2 Loop 3 find descent direction p by solving (4.10); 4 compute step length α for descent direction -p using a line search algorithm; = x + αp; return: x Note, that in step 4 of the Algorithm 6 either backtracking or secant line search can be used. Let us provide numerical simulations to compare these two possibilities and to test the robustness property. The results are presented on Figure 4.4 and Figure 4.5, where the meaning of curves are the same as for Figure 4.1. The system (4.10) is solved by means of Cholesky decomposition.Recall, that the main interest is to see the behavior of the algorithm with and without the presence of noise. The thing to pay attention is how the value of metric behaves when there is no possibility to improve an approximation due to the presence of noise in the model.

Figure 4 . 4 :

 44 Figure 4.4: (a) Metric dist norm between approximation x (k) computed by Algorithm 6 with line search Algorithm 4 and a solution for n = 8, m = 32, and σ = 0, and (b) for σ = 0.1.

Figure 4 . 5 :

 45 Figure 4.5: (a) Metric dist norm between approximation x (k) computed by Algorithm 6 with line search (4.4) and a solution for n = 8, m = 32, and σ = 0, and (b) for σ = 0.1.

Figure 4 . 6 :

 46 Figure 4.6: (a) Profile for solved problems by selected time for 100 problems of size n = 8, m = 32, and σ = 0, and (b) for σ = 0.1. GN b and GN s means Gauss-Newton with backtracking and secant line search respectively.

Proof.m 1 2 m

 2 Let α j denotes the j-th component of Ax, and b j denotes the j-th component of b. Then (4.16) writes as min y∈C j=1 |α j -y j | 2 s.t. |y| = b. (4.17)

Figure 4 . 7 ,Figure 4 . 7 :

 4747 Figure 4.7: (a) Metric dist norm between approximation x (k) computed by Algorithm 7 and a solution for n = 8, m = 32, and σ = 0, and (b) for σ = 0.1.

Figure 4 . 8 :

 48 Figure 4.8: (a) Metric dist norm between approximation x (k) computed by Algorithm 8 and a solution for n = 8, m = 32, and σ = 0, and (b) for σ = 0.1.

Figure 4 . 9 : 1 . 7

 4917 Figure 4.9: (a) Metric dist norm between approximation x (k) computed by Algorithm 9 and a solution for n = 8, m = 32, ρ = 0.5, and σ = 0, and (b) for σ = 0.1.

Figure 4 . 10 :

 410 Figure 4.10: Set of possible projections AA † z for the case of dimension 1 depending on ρ.

Figure 4 . 11 :

 411 Figure 4.11: Set of possible projections AA † z defined by ρ = 1 for the case of dimension 1.

Figure 4 .Figure 4 . 13 :

 4413 Figure 4.12: (a) Metric dist norm between approximation x (k) computed by Algorithm 10 and a solution for n = 8, m = 32, γ = 0, and σ = 0, and (b) for σ = 0.1.

Lemma 1 [ 43 , 1 √ 2

 14312 Lemma 20]. Let a = Re[a] + i Im[a] such that Re[a], Im[a] ∼ N 0, , then for any fixed vector x ∈ C n it holds thatE |a * x| 2 aa * = x 2 I + xx * .

Theorem 8 [ 7 ,

 87 Theorem 3.3] Let x ∈ C n be an arbitrary vector and b = |Ax| ∈ R m + be m samples with m ≥ c 0 • n log n, where c 0 is a sufficiently large numerical constant. Then the computed by Algorithm 12 estimate x (0) normalized to have squared Euclidean norm equal to m

Algorithm 13 :

 13 Gao & Xu initialization [13, Algorithm 1].

  (d) ADMM (Algorithm 10)

Figure 4 . 14 :

 414 Figure 4.14: Time profiles for the phase retrieval algorithms with different initialization strategies: [RI] -random, [WI] -Wirtinger, [GXI] -Gao & Xu.

Figure 4 . 15 :

 415 Figure 4.15: Profiles computed for n = 8, m = 32, σ = 0 and 100 problems in total.

Figure 4 . 16 :

 416 Figure 4.16: Profiles computed for n = 32, m = 128, σ = 0 and 100 problems in total.

Figure 4 . 17 :

 417 Figure 4.17: Profiles computed for n = 128, m = 512, σ = 0 and 100 problems in total.

Figure 4 . 18 :

 418 Figure 4.18: Profiles computed for n = 8, m = 32, σ = 0.1 and 100 problems in total.

Figure 4 . 19 :

 419 Figure 4.19: Profiles computed for n = 32, m = 128, σ = 0.1 and 100 problems in total.

Figure 4 . 20 :

 420 Figure 4.20: Profiles computed for n = 128, m = 512, σ = 0.1 and 100 problems in total.

Figure 4 .

 4 23 and Figure 4.26.

Figure 4 . 21 :

 421 Figure 4.21: Profiles computed for gradient descent algorithm (Algorithm 5) with σ = 0, 100 problems, and Wirtinger initialization (Algorithm 12)

Figure 4 . 22 :

 422 Figure 4.22: Profiles computed for alternating projections algorithm (Algorithm 7) with σ = 0, 100 problems, and Gao & Xu initialization (Algorithm 13)

Figure 4 . 23 :

 423 Figure 4.23: Profiles computed for ADMM algorithm (Algorithm 10) with σ = 0, 100 problems, and random initialization.

Figure 4 . 24 :

 424 Figure 4.24: Profiles computed for gradient descent algorithm (Algorithm 5) with σ = 0.1, 100 problems, and Wirtinger initialization (Algorithm 12)

Figure 4 . 25 :

 425 Figure 4.25: Profiles computed for alternating projections algorithm (Algorithm 7) with σ = 0.1, 100 problems, and Gao & Xu initialization (Algorithm 13)

Figure 4 . 26 :

 426 Figure 4.26: Profiles computed for ADMM algorithm (Algorithm 10) with σ = 0.1, 100 problems, and random initialization.

Algorithm 14 : 2 . 3 .. 2 ) 4 .. 3 ) 5 .

 14232435 Opto-numerical phase correction algorithm [31, Algorithm 1] 1. Initialization: Let x be the known target signal. Let x (0) be the unknown vector of initial phases, such that |x (0) | = | x|. Set k = 0. Magnitude measurement: Let b (k) = |Ax (k) |. Phase retrieval inner algorithm: Compute x (k) ∈ C n as an approximate solution of |Ax| = b (k) .(5Phase modulation: Modify the phases of input signal according tox (k+1) = x (k) exp(i(arg( x) -arg( x (k) ))).(5Return loop: Set k = k + 1 and go to 2.

Figure 5 . 1 :

 51 Figure 5.1: Time profiles for the phase retrieval algorithms to solve phase correction problems with different initialization strategies: [RI] -random, [WI] -Wirtinger, [GXI] -Gao & Xu.

Figure 5 . 2 :

 52 Figure 5.2: Profiles computed for n = 8, m = 32, σ = 0 and 100 problems in total.

Figure 5 . 3 :

 53 Figure 5.3: Profiles computed for n = 32, m = 128, σ = 0 and 100 problems in total.

Figure 5 . 4 :

 54 Figure 5.4: Profiles computed for n = 128, m = 512, σ = 0 and 100 problems in total.

Figure 5 . 5 :

 55 Figure 5.5: Profiles computed for n = 8, m = 32, σ = 0.1 and 100 problems in total.

Figure 5 . 6 :

 56 Figure 5.6: Profiles computed for n = 32, m = 128, σ = 0.1 and 100 problems in total.

Figure 5 . 7 :

 57 Figure 5.7: Profiles computed for n = 128, m = 512, σ = 0.1 and 100 problems in total.

Figure 6 . 1 :

 61 Figure 6.1: Average q norm after K max = 8 corrections for n = 16, m = 64.

Figure 6 . 2 :

 62 Figure 6.2: Average q norm after K max = 1 corrections for different n.

Figure 6 . 3 :

 63 Figure 6.3: Average q norm after different K max corrections for n = 16, m = 64.

Conjecture 1

 1 Let F(x) : C n → R m + be such that F(x) = |Ax|. Let us consider an abstract inverse function F -1 (x) : R m + → C n which solves (5.2) exactly. Let NN(x) be a function of the form (6.1) trained with the Algorithm 15 for a fixed target x ∈ C n . Let B ε ( x) ⊂ C n be a ball centered at x with radius ε > 0. Then q norm x, e i T F N N (b) ≈ 0 where b = F(x) for any x ∈ B ε ( x).

Remark 2

 2 It was observed numerically that if | x| is not a constant then the sufficiently small values of q norm can not be reached during training. (see Figure6.4).

Figure 6 . 4 :

 64 Figure 6.4: Average q norm after different K max corrections for n = 16, m = 64, where red dashed line is a required threshold q norm = 0.04.

Figure 6 . 5 :

 65 Figure 6.5: Average q norm after K max = 8 corrections for n = 16, m = 64 during training, where red dashed line is a required threshold q norm = 0.04.

Observation 9 2 Observation 10 2 Proposition 14

 9210214 Let λ > 0 and x ∈ R m , then TFNN(λx) = TFNN(x).Proof. Let W ∈ C n×m , then by linearity arg(W (λx)) = arg(λW x) = arg(W x) 1 . Let λ > 0, A ∈ C m×n , and x ∈ C n , then |A(λx)| = λ|Ax|. Proof. |A(λx)| = |λAx| = |λ||Ax| = λ|Ax|. Let x ∈ C n be such that |x| = λ for λ > 0. Let A ∈ C m×n . Then TFNN |Ax| = TFNN |Ae i arg(x) | .Proof. Using the Observation 10 we have that TFNN |A(λe i arg(x) )| = TFNN λ|Ae i arg(x) | , and then with the Observation 9 we get TFNN λ|Ae i arg(x) | = TFNN |Ae i arg(x) | which shows the invariance to a constant amplitude.2

Observation 11 2 Proposition 15

 11215 Let A ∈ C m×n be a transmission matrix. Then, |(A e -i arg(a •1 ) )x| = |Ax| for any x ∈ C n .Proof. Let us consider the j-th row of A and show that |(a j e -i arg(a j1 ) ) x| = |a j x|.|(a j e -i arg(a j1 ) ) x| = n k=1 a jk e -i arg(a j1 ) x k = e -i arg(a j1 ) n k=1 a jk x k = n k=1 a jk x k = |a j x|.Applying this result to each row we finish the proof. Let A ∈ C m×n be a transmission matrix and let A its normalized form (see Definition 5 in Appendix A). Then TFNN(| Ax|) = TFNN(|Ax|). Proof. Applying the Observation 9 and the Observation 11 we obtain TFNN(| Ax|) = TFNN 1 A ∞ (A e -i arg(a :1 ) )x = TFNN (A e -i arg(a :1 ) )x = TFNN(|Ax|).

Figure 6 . 6 :

 66 Figure 6.6: (a) Heat maps of the minimal achievable mean q norm in grey scale and (b) its required relative training time. The relative time on (b) is computed by dividing a learning time in seconds for each n and m/n by the minimal time to obtain GPU invariant information. The minimal time required by the GPU used in for this experiments was 1.3 s.

Figure 6 . 7 :

 67 Figure 6.7: Mean q norm behavior during corrections.

Figure 6 . 8 :

 68 Figure 6.8: Average q norm after K max = 8 corrections for n = 6, m = 64.

Figure 7 . 1 :

 71 Figure 7.1: Profiles computed for n = 16, m = 64, 100 random starting lasers state.

Figure 7 . 2 :

 72 Figure 7.2: q norm experimental traces computed for n = 16, m = 64, 100 random starting lasers state.

Figure 7 . 3 :Figure 7 . 4 :

 7374 Figure 7.3: Profiles computed for n = 16, m = 192, 100 random starting lasers state.

Figure 7 . 5 :

 75 Figure 7.5: Experimental setup picture.

Algorithm 17 :

 17 Transmission matrix computation. Input: Signals X ∈ C N ×n and measurements B ∈ R N ×m + , tolerance ε > 0, maximal number of restarts R max ∈ N\{0}. Output: An approximation of the transmission matrix A ∈ C m×n . 1 Initialize each row of matrix A (0) ∈ C m×n by solving m phase retrieval problems for X and B with the starting point defined by Gao & Xu initialization [13]; 2 Set k = 0; 3 for k < R max do 4

5 if J is empty then 6 break 7 Find

 567 an approximate solution a j ∈ C n of equations |Xa| = b •j for all j ∈ J by means of Algorithm 7, where the starting points are generated randomly; 8 Build the matrix A (k+1) where its rows a k = k + 1;9 return A = A (k) .

  Figure A.2 it can be clearly seen that unbiasing technique improves the quality

Figure A. 2 :

 2 Figure A.2: Comparison of (A.3) for biased and unbiased models using experimental data X and B for n = 16, m = 64 and N = 1000.

Figure A. 3 :

 3 Figure A.3: Traces of q norm during experimental phase corrections performed by the alternating projections algorithm.

Definition 5

 5 Let A ∈ C m×n be a transmission matrix. Then the normalized matrix A computes asA = 1 A(:) ∞ A e -i arg(a •1 ) , (A.4)where A(:) ∞ = max i∈{1,...,m},j∈{1,...,n} |a ij | and a •1 = (a 11 , a 21 , . . . , a m1 ) is the first column of A.

  Let z = (z 1 , . . . , z n ) ∈ C n be an n-dimensional column vector, where z j = x j + iy j for j ∈ {1, . . . , n}. Let f : C n → R be a real-valued function of complex variable z ∈ C n Then, Wirtinger gradient writes as

									∂z	(Conjugation rule)
		df =	∂f ∂z	dz +	∂f ∂ z dz	(Differential rule)
		∂(h • g) ∂z	=	∂h ∂g	∂g ∂z	+	∂h ∂ḡ	∂ḡ ∂z	(Chain rule)
		∂(h • g) ∂ z	=	∂h ∂g	∂g ∂ z +	∂h ∂ḡ	∂ḡ ∂ z	(Chain rule)
	Observation 2 (Wirtinger gradient [16]) ∇f (z) = 2	∂f ∂z	,	∂f ∂ z	*	,
	where * operation means consequent conjugation and transposition, and ∂f ∂z = ∂f ∂z 1 , . . . , ∂f ∂zn	,
	∂f ∂ z = ∂f ∂ z1 , . . . , ∂f ∂ zn	for ∂f ∂z j = 1 2	∂f ∂x j -i ∂f ∂y j , ∂f ∂ zj = 1 2	∂f ∂x j + i ∂f ∂y j .
	Proof. Let us consider a coordinate transformation (2.6), which we augment to n-dimensional
	case as

  Let the total number of iterations is set to 100. Let the tolerance value is set to 0.001. Let the distance metric is set to dist norm . Let x (k) ∈ C n be the k-th approximation of a solution computed by some algorithm, and let x ∈ C n be a solution. Then, the problem is considered to be solved if for all k > 90 we have that dist norm (x (k) , x) < 0.001.

	2.5. PROFILING METHOD
	Example 2 (Solved problem) Definition 2 (Time measurement) The time of one iteration of an algorithm is measured
	by counting the number of performed basic operations and multiplying them by their correspond-
	ing time of evaluation.
	Definition 3 (Time of basic operations) The list of basic operations and their correspond-
	ing time of evaluation, which is CPU dependent and fixed in this work, is presented in Ta-
	ble 2.1. These values were obtained by measuring the time of performing the selected operations
	Operation Time in seconds
	+	0.1089426 • 10 -8
	-	0.0834372 • 10 -8
	×	0.1665280 • 10 -8
	/	0.3297135 • 10 -8
	| • | √ •	0.2500649 • 10 -8 0.1788236 • 10 -8
	sin(•)	1.4113816 • 10 -8
	cos(•)	1.1528648 • 10 -8
	tan(•)	2.5544987 • 10 -8
	arctan(•) 2.2179606 • 10 -8
	exp(•)	1.3695341 • 10 -8

Table 2 .

 2 1: Time to perform basic operations.

  and the unknown is the couple (y, z) ∈ C n × C m . Obviously, x is a solution of (1.2) if and only if (Ax, Ax) is a solution of (4.20).The system of equations (4.20) can be reformulated as an artificial minimization problem of a zero function subject to set and equality constraints:min (y,z)∈C m ×C m 0 s.t.Let us consider the minimization problem (4.21). Then the Altenating Directions Method of Multipliers update formulas are

y -z = 0, y ∈ range(A), z ∈ M b (4.21)

Let us define the update formulas for the Altenating Directions Method of Multipliers which is applied to (4.21).

Proposition 8

  The best initialization strategy for ADMM is Gao & Xu initialization[START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF]. Note, that it is different from what we observe on Figure4.14d. The reason comes from different metric. It reveals that GXI is better to use if q norm is considered as a distance measure, and RI if dist norm .

	2. GXI and WI initialization strategies are equivalent for gradient descent algorithm (Fig-
	ure 5.1a).
	3. The best initialization strategy for alternating projections algorithm is Gao & Xu initial-
	ization [13].
	4.

1.2. PHASE RETRIEVAL PROBLEM

In[START_REF] Saucourt | Random scattering and alternating projection optimization for active phase control of a laser beam array[END_REF] the phasing quality is defined as 1 -q(x, y)

.

3.7. QUASI-REINFORCEMENT LEARNING

return: x = x (k+1)

The same results holds with the intialization from Algorithm 12 because m i=1 a i n ≈ mn with a standard deviation of about the square root of this quantity.

Geometrically, it can be considered that the angle of vector does not change if we multiply it by the positive constant.

Remerciements

Target fixed neural network

In this section we present the neural network based algorithm to solve the phase retrieval problem in step 3 of the Algorithm 14 approximately for a fixed target x.

Let TFNN(x) : R m → [-π, π] n be a neural network model such that TFNN(x) = arg(W x), (

where W ∈ C n×m is a complex matrix of trainable parameters. The Algorithm 15 was developed to find W such that TFNN can compute phase corrections. Algorithm 15: QRL algorithm for TFNN 1. Initialization: Let N ∈ N be the size of the batch. Let K max ∈ N be a number of corrections. Let I max ∈ N be a number of iterations. Let x ∈ C 1×n be the known target signal. Let W ∈ C n×m be a random matrix. Set l = 0.

2. Data generation: Generate a batch of initial phases X (0) ∈ C N ×n , such that |x (0)

(Step 3b in Algorithm 2 where ∆w

applying the Adam update rule [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with parameters α = 0.001, β 1 = 0.9, β 2 = 0.999. (Step 3c in Algorithm 2)

7. Phase modulation: Modify the phases of input signal according to The idea is to plug the model in the simulated phase correction loop in order to train the network to converge to target from any starting signal x (0) . This idea can be formulated in terms of the quasi-reinforcement learning algorithm (Algorithm 2) where the agent Agent is the model (6.1), the action Ât is a phase correction vector ϕ (k) = arg( x) -arg( x (k) ), the state change appears. In this section, we present the neural network to use in the phase correction loop where a target changes with a high frequency. The matrix W in (6.1) implicitly depends on x, which means that there exists a mapping TANN : [-π, π] n → C n×m . Let us construct this mapping in the following form TANN(x) = Reshape(U e ix ), (6.4) where U ∈ C nm×n is a complex matrix of trainable parameters and Reshape : C mn → C n×m is a matrix from vector building operator. The approach to learn (6.4) is presented in the Algorithm 16.

Algorithm 16: QRL algorithm for TANN 1. Initialization: Let M ∈ N be the size of the batch for signals and N ∈ N for targets. Let K max ∈ N be a number of corrections. Let I max ∈ N be a number of iterations. Let U ∈ C nm×n be some random matrix. Set l = 0.

2. Data generation: Generate a batch of initial targets X ∈ C N ×n and initial phases

Parameters update: Update U applying the Adam update rule [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with parameters α = β 1 = β 2 = 0.1.

8.

Phase modulation: Modify the phases of input signal according to

10. Return iteration loop:

Let us explain each step of the Algorithm 16 in details. To simplify the notations for explanations of steps 2-8 let us fix k and remove it from the notations, and consider that j ∈ {1, . . . , N }, p ∈ {1, . . . , M } where it is not explicitly specified.

In step 1 the parameters M , N , K max , I max and U are initialized. The information about K max and I max is presented in parts in the explanation of the Algorithm 15. The parameters N and M represent the size of training batch. In contrast to the Algorithm 15, it is necessary to generate not only initial signals to correct, but also targets. The parameter N represents was recorded and visualized as a heat map in Figure 6.9a, with the corresponding relative training time shown in Figure 6.9b. The minimal achievable q norm is obtained by solving 1000 phase correction problems with different targets for each combination of n and m/n. The red line in Figure 6.9a reveals the dependency between n and m/n to obtain q norm = 0.04 and is defined as f (n) = n 2 + 1. This gave us information about the minimal number of measurements needed to obtain q norm ≤ 0.04, which was m = n 2 2 + n. An average q norm during corrections is visualized in Figure 6.10, which gives the information about the number of corrections required to achieve a target.

Appendix A Transmission matrix computational details

In this section, we give additional details on how a transmission matrix A ∈ C m×n can be computed to build a mathematical model (1.2), which then is used for phase corrections in Algorithm 14. As it is explained in the Introduction chapter, we obtain matrix A by solving the following problem

where the sets of signals X ∈ C N ×n and measurements B ∈ R N ×m + are obtained experimentally. This problem can be split into m phase retrieval problems that can be solved by one of the methods presented in this work. In addition, there are no time constraints. For this goal, we select the alternating projections algorithm since it does not require any additional parameters that sensitive to the level of noise like γ in ADMM. Recall, that even with Gao & Xu initialization strategy [START_REF] Gao | Phaseless recovery using the gauss-newton method[END_REF] there is no guarantee that the alternating projections algorithm will find a global solution. Since this part is critical and directly affects the quality of the model, algorithm with restarts is considered, where the main idea is to change a starting point until a global solution is found. However, it is not obvious how to determine the global solution. In practice, there is no A such that |XA | = B nearly exactly because of the presence of noise. The difficulty is that the level of noise is not known in advance and thus the choice of a threshold is not obvious. That is why, we propose to use the following statistical rule to answer this question. Let A be an approximate solution of |XA | = B. Then the j-th row a j of A is a global solution of

Here, µ and s are sample mean and standard deviation of all m relative norms. μ is an upper bound of a trust region for true value of µ with a critical value z α of the standard normal distribution for the significance level α. We set α = 0.05 which means that the true value of µ is and modeled measurements as

where a typical value is ≤ 0. dependence. For this goal, we present the Algorithm 18 which can compute the shift value, which must be applied to all measurements before transmission matrix computation.