Real-time Adaptive Optics with pyramid wavefront sensors: Accurate wavefront reconstruction using iterative methods

In this paper, we address the inverse problem of fast, stable, and high-quality wavefront reconstruction from pyramid wavefront sensor data for Adaptive Optics systems on Extremely Large Telescopes. For solving the indicated problem we apply well-known iterative mathematical algorithms, namely conjugate gradient, steepest descent, Landweber, Landweber-Kaczmarz and steepest descent-Kaczmarz iteration based on theoretical studies of the pyramid wavefront sensor. We compare the performance (in terms of correction quality and speed) of these algorithms in end-to-end numerical simulations of a closed adaptive loop. The comparison is performed in the context of a high-order SCAO system for METIS, one of the first-light instruments currently under design for the Extremely Large Telescope. We show that, though being iterative, the analyzed algorithms, when applied in the studied context, can be implemented in a very efficient manner, which reduces the related computational effort significantly. We demonstrate that the suggested analytically developed approaches involving iterative algorithms provide comparable quality to standard matrix-vector-multiplication methods while being computationally cheaper.


Introduction
We consider the problem of fast, stable and highly accurate wavefront correction for large-scale real-time closed loop Adaptive Optics (AO) systems on Extremely Large Telescopes (ELTs). More specifically, we focus on reconstructing the wavefront Φ from pyramid wavefront sensor data s = [s x , s y ] with a pre-defined accuracy and within the required time slot to fit the real-time setting. The available measurements s are related to the wavefront Φ by the non-linear integral operator P = [P x , P y ] representing the pyramid sensor model. In part I of this paper [36], we extensively studied the physical and mathematical forward models of the pyramid wavefront sensor and the underlying operators. The major aim of this paper is to apply the theory of part I [36] in order to solve the inverse problem of reconstruction the unknown wavefront from the given sensor data. Several approximations of the full Fourier optics based, non-linear pyramid sensor model derived in [36] allow for a fast numerical implementation of the wavefront sensor (WFS) operators, which suggests to apply iterative algorithms for solving the reconstruction problem.
We study and compare the performance (in terms of correction quality and speed) of well-known mathematical algorithms for solving inverse problems, namely conjugate gradient, steepest descent, Landweber, Landweber-Kaczmarz, and steepest descent-Kaczmarz iteration. The comparison is performed within the context of two specific high-order AO systems which are currently under design for the future ELT instruments. One of the systems is the SCAO (Single Conjugate Adaptive Optics) module for the METIS instrument [6] on a 39 m telescope equipped with a 74 × 74 pyramid wavefront sensor (PWFS) sensing in the near-infrared K-band (at the wavelength λ = 2200 nm). The SCAO system is supposed to control ∼ 4000 mirror actuators at frequencies of 500 − 1000 Hz. The second system of interest is the eXtreme Adaptive Optics (XAO) module for the EPICS instrument [44] having a 200 × 200 PWFS as a core component. For the XAO system a huge amount of ∼ 30000 mirror actuators have to be controlled at a frequency of 3 kHz, which is a challenge for the real-time control as the algorithm has to produce the reconstruction within a time of less than 0.33 ms. Both considered systems are expected to provide perfect correction quality resulting in an unprecedented image contrast required for their scientific aims.
Additionally, we compare the performance of our methods to the so far standard wavefront reconstruction algorithms used in these days. These are based on matrix-vector-multiplication (MVM) and invert the most exact Fourier optics model. However, the main drawback of any interaction-matrix-based method is their computational effort, which is demanding for the planned large-scale real-time AO systems. The computational complexity required for setting up the command matrix scales as O(N 3 ) with N being the number of controlled actuators, and the application of this command matrix on the sensor data scales as O(N 2 ) [17]. The indicated limitation of MVM methods makes their application on large-scale AO systems, such as the XAO system on the ELT, hardly feasible even on the hardware expected at the time of the telescope launch in around 2024. The steadily growing mirror sizes of future telescope systems imply an immense grow in the computational load of existing algorithms which means that the numerical effort has to be kept in mind when developing new algorithms. We demonstrate in this paper that all the proposed iterative algorithms provide the required high-quality and stable wavefront correction along with the heavily reduced computational complexity compared to standard matrix-vector-multiplication methods. The suggested matrix-free approaches make highly accurate real-time wavefront reconstruction feasible even for the XAO system.
Since the pyramid wavefront sensor was first introduced in astronomical AO by Ragazzoni [60], it has been extensively studied in optical test benches [21,58,61,78] and acknowledged to possess several precious characteristics distinguishing it from other types of wavefront sensors. Among those are the increased sensitivity, adjustable linearity range, and pupil sampling, as well as its ability to sense the segmented piston modes [21,58], which is gaining a special importance in the era of ELTs inevitably having segmented mirrors. Due to the named advantages, the PWFS is nowadays integrated as baseline on several telescope instruments under design. Among those we can name MICADO [12], HARMONI [27,55], METIS [6], EPICS [44], and ATLAS [26] on the ELT, SCAO and LTAO systems on the Giant Magellan Telescope (GMT) [76,22], as well as the NFIRAOS [53,77] and PFI [52] instruments on the Thirty Meter Telescope (TMT). Furthermore, the application area of the PWFS is not limited to Adaptive Optics in Astronomy. Apart from astronomical observations, the pyramid sensor is also utilized in AO in ophthalmology [10,15,16,40] and microscopy [39,41]. Hence, the problem of improving the quality of wavefront reconstruction approaches for this type of sensor by using sophisticated mathematical methods has never been more interesting, challenging, and important.
In Section 2 we describe algorithms existing for wavefront reconstruction from pyramid sensor data. We sketch their advantages, drawbacks and limitations. We proceed with recalling the theoretical principles of wavefront sensing using the pyramid sensor in Section 3 and mention details on the discretization of the sensor in Section 4. Afterwards, we describe several iterative algorithms, namely the conjugate gradient method for the normal equation and the steepest descent method, Landweber iteration as well as Kaczmarz type algorithms in Section 5. This Section contains details on the numerical implementation of the involved operators as well. Section 6 presents the performance of the proposed algorithms and a comparison with respect to the achieved reconstruction quality. Finally, in Section 7 we evaluate and compare the computational complexities of the analyzed approaches. Both, the reconstruction quality and the speed of the algorithms are additionally compared versus those of an MVM approach.

Existing algorithms
For the PWFS, several types of reconstruction approaches have been considered so far [33,38]. Frequently used are interaction-matrix-based methods which were already mentioned in Section 1. A more computationally efficient algorithm, the Fourier Transform Reconstructor (FTR), with a complexity of O(N log N ), was suggested in [59]. The method was developed for a simplified geometrical optics based model of the roof WFS, which assumes a large amount of modulation applied to the sensor. Under these assumptions the sensor data are modeled as the derivative of the phase, similar to the Shack-Hartmann (SH) sensor. In the reported algorithms, the authors applied SH Fourier domain filters for wavefront reconstruction. A slightly worse performance of the SH-based FTR compared to MVM was demonstrated for an 8 m telescope. However, since the geometrical model is valid only for large modulations, this method requires further research on a modal optimization, which becomes especially important in the case of ELTs. The idea of applying inverse Fourier domain filters for wavefront reconstruction from pyramid sensor was further developed in [72,75], where the authors reported on two methods -Convolution with the Linearized Inverse Filter (CLIF) with the complexity O(N 3/2 ) and the Pyramid Fourier Transform Reconstructor (PFTR) with O(N log N ) complexity. Both algorithms use a precise Fourier optics forward model allowing the derivation of more exact Fourier filters connecting the incoming phase to the sensor measurements. The fastest available reconstruction algorithm is the Preprocessed Cumulative Reonstructor with Domain decomposition (P-CuReD) [73,74]. Based on an analytical relation in the Fourier domain, pyramid data are transformed into SH-like data and subsequently inverted using CuReD [65,66] which was originally developed for SH sensors and has already been tested on-sky [3,4]. This algorithm provides in numerical simulations the same or even better quality results than the standard interaction-matrix-based approaches and scales with a complexity of O(N ). Further algorithms are based on the fact that the wavefront reconstruction from the full pyramid sensor model can be simplified to an inversion of the finite Hilbert transform. The FHTR (Finite Hilbert Transform Reconstructor) algorithm proposed in [71] uses the direct inversion formula of the finite Hilbert transform and a second method called SVTR (Singular Value Type Reconstructor) [35] is based on an analytical singular value expansion of the finite Hilbert transform operator. Both algorithms have the complexity O(N 3/2 ).

Pyramid/Roof wavefront sensors
In the following, we recall the theory discussed in part I of the paper [36]. Since all algorithms are based on an approximation of the full pyramid sensor, or more precisely on a linearization of the roof sensor, we will only focus on the properties of the latter and mention corresponding adjoint operators which are needed for the application of the proposed iterative methods. For a precise analysis of the full pyramid sensor model, details about the linearization procedure and proofs we refer the reader to part I of the paper [36]. The analytical models of the non-and modulated pyramid wavefront sensors are complex and therefore difficult to invert directly. However, some assumptions suggested by the physical setting of the model itself, allow to simplify the non-linear Fourier optics based model of the pyramid sensor P = [P x , P y ] (cf [36] for the definition) by substituting the pyramidal prism by two orthogonally placed two-sided roof prisms [7,57,78]. Due to the physical decoupling of the roof prisms and their orthogonality, the two signal sets s = [s x , s y ] provided by the pyramid sensor as become independent and contain information about the incoming phase Φ only in x-and y-direction respectively. Because of symmetry, we only consider the roof sensor operator R = [R x , R y ] in x-direction, i.e., R x . By interchanging x and y all assertions are derived for R y accordingly. We describe the annular telescope aperture mask by Ω = Ω y ×Ω x ⊆ [−D/2, D/2] 2 . Single lines (intervals) of the annular aperture are represented by Ω x = [a x , b x ] and Ω y = [a y , b y ], with a x < b x , a y < b y being the borders of the pupil for fixed x and y correspondingly. with R {n,c,l},lin

Linearized roof WFS
indicating the linearized roof sensor operators and The modulation parameter is given by α λ = 2πα λ and J 0 denotes the zero-order Bessel function of the first kind.
The linearized roof sensor operators R {n,c,l},lin offer a further possibility for simplification of the model due to the splitting for the integral operators L {n,c,l} where p.v. denotes the Cauchy principal value. Note that L n x is the finite Hilbert transform operator. Dropping the second term in (3) leads to the inverse problem for pyramid sensor data s x .

Adjoint operators
All iterative methods for wavefront reconstruction from pyramid sensor data proposed below require the application of adjoint operators. As mentioned, e.g., in [36], it holds that Φ ∈ H 11/6 , i.e., the underlying operators are defined as L {n,c,l},lin Proof. See [36].

The discrete pyramid wavefront sensor
The full continuous measurements s x (x, y) and s y (x, y) of the pyramid wavefront sensor are not available in practice. For the description of the discrete pyramid sensor we perform a division of the continuous two dimensional process into finitely many equispaced regions called subapertures. The data are then assumed to be averaged over every subaperture which corresponds to the finite sampling of the pyramid sensor. Note that in reality, the subaperture grid is predefined by the sensor's physics. Following the approach in [78], we examine the sensor data as functions evaluated in the (discrete) middle points of the WFS subapertures. In the two dimensional case we consider quadratic subapertures of size d × d with d = D n , where D represents the telescope diameter, i.e., the primary mirror size, and n the number of subapertures in one direction.
Note that all considerations are valid for measurements both in x-direction s x (x, y) and y-direction s y (x, y), as well as for non-modulated, circularly, and linearly modulated data. Thus, we consider general measurements identified by s(x, y). Discretizing s(x, y) delivers n 2 data values s jk with j, k = 1, . . . , n.
For the following, we use the Dirac comb III d defined as where δ denotes the delta distribution.
The continuous signal is captured by the wavefront sensor as follows: First, the average of the measurements over one subaperture is calculated. This is represented as a convolution of the continuous data s(x, y) with a characteristic function The discretization is carried out as a application of the Dirac comb III d assuming that the measurements s fulfill the necessary conditions on applying the distribution δ. Herewith, we assign a discrete set of measurements s centered on the subapertures from floating average valuess(x, y) to the discrete set of subaperture middle points {(jd, kd) : j, k ∈ Z}.
Finally, we restrict the number of measurements to the size of the region captured by the sensor. For several telescope systems the pupil Ω is annular instead of circular since a shade created by the secondary mirror prevents measurements on these areas. Thus, the light in the area of the central obstruction possibly does not produce reliable measurements [20,37,56,69,70]. Hence, the last step is realized as a multiplication with a second characteristic function for j, k = 1, . . . , n. To be precise, we usually consider less than n 2 measurements due to the annular shape of the aperture and ignore subapertures which are too less illuminated in order to produce reliable data. However, we will not specifically mention this fact throughout the paper.

Iterative wavefront reconstruction methods
In this Section, we adapt well-known mathematical algorithms, namely the conjugate gradient method for the normal equation, the steepest descent algorithm, and Landweber iteration as well as modifications of these methods coupled with a Kaczmarz strategy to the problem of wavefront reconstruction from pyramid sensor data. By solving the WFS equations in the operator setting instead of forming matrices we minimize the computational complexity of the proposed methods. For wavefront reconstruction we solve the two integral equations representing a linearization of the roof WFS according to (2) and s = [s x , s y ] pyramid sensor measurements. As a further simplification, we consider the inverse problem For simplicity of notation we use Q : in the following since the basic idea is the same for all types of modulation. Moreover, we only concentrate on solving the inverse problem (5)-(6), but mention that solutions of (1) and (7) can be calculated accordingly.
To specify the representation of the incoming phase Φ and the measurements s we denote the number of subapertures by n. There are various possible representations for the phase and the measurements, e.g., Zernike polynomials or bilinear spline functions. We choose a representation that guarantees maximum computational efficiency, and thus assume that the incoming phase and the measurements are piecewise constant on the subapertures, i.e., where (φ i ) 1≤i≤n , (s x,i ) 1≤i≤n denote basis coefficients and Ω y i = x y i−1 , x y i the i-th subaperture of a row for fixed y. As the wavefront sensor provides two measurements (one in x-and one in y-direction), for every single subaperture, the suggestion of representing the measurements via piecewise constant functions describing the subaperture grid is reasonable. We calculate the involved operators as, e.g., The functions α {n,c,l} i (x, y) are computed offline and do not influence the computational speed of the proposed methods. The implementation of all involved operators is performed analogously when choosing the basis representation (8). Now, we focus on concrete wavefront reconstruction algorithms for pyramid sensor data using iterative methods. In particular, we consider the conjugate gradient algorithm for the normal equation (CGNE), the steepest descent (SD), and Landweber iteration. Since the pyramid sensor provides two measurements s x and s y , the above named approaches deliver two solutions Φ = [Φ x , Φ y ], one in x-and one in ydirection. The final reconstruction Φ rec is then computed as the average of the two temporary solutions. An alternative combination of the two measurements s x and s y using Kaczmarz loops is investigated in Section 5.4. Note that the theory of the presented algorithms is mainly based on [19,51]. The considered norms are the L 2 -norms.

CGNE approach
The conjugate gradient (CG) method is one of the most powerful algorithms for solving self-adjoint, positive (semi-)definite linear equations [5,19,29,30,31,43]. For solving the wavefront reconstruction problem we apply the conjugate gradient method to the normal equation [19] by requiring the fewest iterations among all semiiterative methods.
The CGNE method (Algorithm 1) applied to the inverse problem of wavefront reconstruction from pyramid wavefront sensor data is described by:

Steepest descent approach
For solving the system (5)-(6) we are additionally interested in the method of steepest descent where we consider different choices of the step sizes in the iterative process. For pyramid sensors, we use the SD method (Algorithm 2) applied to the least-squares functional The method of steepest descent was originally introduced by Cauchy [8] as one of the most basic procedures to minimize a differentiable functional. A popular step size is determined by an exact line search in the direction of the negative gradient. Alternative choices of the step size have already been considered, e.g., in [67,68] and will be discussed below for the problem of wavefront reconstruction from pyramid data using the SD method. The gradient of the classical least-squares functional is given by and the resulting algorithm reads as:

Step size choices and convergence
The speed of convergence of the gradient iteration depends highly on the choice of the step size τ i . We consider the classical steepest descent (line search) step size that is defined by This means that an exact line search is performed in the direction of steepest descent which corresponds to the direction of the negative gradient. For the least-squares functional (10) and corresponding derivative (11) the steepest descent step size with d i = −Q * (QΦ − s) reads as and results in the so called Cauchy method.
If we minimize the gradient norm along the search direction, we obtain another line search method for finite dimensions, namely the method of minimal gradient (MG) [14] given by Because of zigzagging between consecutive steps the SD method suffers from slow convergence in some cases. To overcome these effects, a fast and efficient alternative step size choice was introduced by Barzilai and Borwein (BB) in [2]. The BB technique is motivated by quasi-Newton methods and derived from a two-point approximation to the secant equation. There exist two versions of the BB method which are defined by Plugging in the calculations corresponding to Algorithm 2 we obtain since the involved operators are linear. Therefore, the BB step sizes are rewritten as i.e., The idea is to use additional information of the previous iteration to compute the step size for the current iteration. Once again with the Cauchy-Schwarz inequality we obtain τ BB1 i ≤ τ BB2 i . Generally, while the SD and MG method decrease monotonically, the BB step size choices are non-monotone as the error behaves non-monotonously, i.e., ||Φ − Φ i+1 || ≤ ||Φ − Φ i || for the true solution Φ is not fulfilled for every iteration i. Nevertheless, the BB method converges to a solution of (10) as found in [63].
The Cauchy-Barzilai-Borwein (CBB) step size is based on the idea to use the SD and BB step size alternating. The method, which was introduced in 2003 and is also called alternate step size (AS) gradient method, aims at reducing the zigzag-effect of the Cauchy method [13], and therefore leads to a faster convergence. The promising alternative to the BB method reads as for i even.
Due to τ BB2 i = τ SD i−1 , we use the same step size twice in two consecutive iterations. An alternate version of the MG method called alternate minimization (AM) gradient method was proposed in [14] having an SD iteration for every second step. Generally, the SD method becomes faster when one non-monotone (e.g., BB) step is made even after several SD steps [79]. In addition, a variety of step size choices have been introduced using combinations or shortened step size versions of the above mentioned options.
To reduce the computational effort of the algorithms, one can use a fixed step size.
Step sizes which are tuned heuristically depend mainly on the size of the telescope and the resolution of the wavefront sensor (discretization). For a fixed step size τ i = β the steepest descent algorithm for the least-squares functional (10) reduces to the standard Landweber iteration.
The Landweber iteration modified for wavefront reconstruction based on pyramid sensor measurements (Algorithm 3) reads as: Besides the above discussed methods for pyramid sensors, there already exist several algorithms providing two reconstructions, one from data s x and one from data s y [35,71,75]. Since in the reconstructions obtained by averaging the two solutions we experienced distinct horizontal and vertical artifacts, the aim is to combine the reconstructions already during the iteration steps. For this reason we investigate Kaczmarz methods in which the two data sets are used alternating.

Kaczmarz methods for wavefront reconstruction from pyramid sensor data
For the reconstruction of the incoming wavefront, the pyramid sensor provides two data sets s x and s y . If the reconstruction were based on the full pyramid model, the incoming phase Φ either could be reconstructed solely from measurements s x or solely from s y because the null space of the operators consists only of the global piston mode, which anyway does not influence the imaging quality. However, in case the reconstruction algorithms utilize the roof sensor model, both data need to be used due to different null spaces of the single operators. Altogether, there are several facts that support the usage of both data sets. On the one hand, we expect better reconstruction quality in case we utilize more information. This argument is additionally strengthened by the presence of noise in the sensor measuring process. On the other hand, deeper investigations of the underlying operators in x-and y-direction show that they have different null spaces, i.e., depending on the underlying model of the reconstructors there exist modes that cannot be reconstructed. For instance, pyramid and roof wavefront sensors are not able to detect a constant added to the incoming phase Φ. This undetectable constant, called piston mode (mode of order 0), has no influence on the measurements s. In order to characterize effects that are invisible in sensor data we discuss selected wavefront modes (of order 0 and 1) which are elements of the null space of the roof wavefront sensor operators, i.e., phase elements that deliver measurements equal to zero. For the following investigations, we will consider the mathematical forward model of the linearized roof sensor R lin and analyze the null spaces of the corresponding operators described by We study the response of the linearized roof sensor to a global piston mode shown in Figure 1 left. Hence, where c ∈ R is a constant. Furthermore, we analyze how the sensor responses to modes of order 1 called Proof. Global phase piston modes c · X Ω are in the null space of the roof sensor operators because of and R {n,c,l},lin y (x, y) respectively.
Then, we consider Altogether, we obtain that tip is in the null space of R {n,c,l},lin y and tilt in the null space of R {n,c,l},lin x . Figure 1: The figures indicate piston, tip and tilt mode (from left to right).
If we further simplify the roof sensor model by excluding the second term in (2) and consider the operators L {n,c,l} defined in (4), we note that these operators are injective, i.e., N L {n,c.l} = {0}. This assertion follows from the injectivity of the finite Hilbert transform shown in [18].
Note that the results of Proposition 2 are directly transferred to the full non-linear roof sensor operators R {n,c,l} discussed in part I of the paper [36]. As soon as we consider functions including various powers of x or y in Φ tip/tilt , the corresponding measurements in x-or y-direction are equal to zero as well. The fact that R {n,c,l},lin x and R {n,c,l},lin y have different null spaces intensifies the requirement of an appropriate combination of the two data sets s x and s y for reconstruction methods which are based on the roof sensor model.
One idea to appropriately combine both data sets is to reconstruct independently in both directions and average the reconstructions at the end as already considered for the CGNE, SD, and Landweber iteration above. However, it is not guaranteed that the final (averaged) solution Φ rec fulfills both equations (5)-(6). Another possibility is two consider Q x , Q y as one single operator and a third one is to use a Kaczmarz strategy [42,54] which is computationally cheaper and for which it is guaranteed that the equations (5)-(6) are fulfilled for the final solution. Kaczmarz methods, in general, have been developed for solving linear systems of equations. We have decided to implement Kaczmarz strategies for the pyramid sensor in combination with several of the above discussed algorithms.

Landweber-Kaczmarz approach
In practice, the Landweber algorithm is used because it is simple and each iteration is cheap. Though, the process usually requires a high number of iterations. Anyway, we do not experience slow convergence for reconstruction from pyramid data due to a close similarity between adjoint and inverse operators as investigated in [35] for the non-modulated sensor, i.e., for the finite Hilbert transform operator. When using proper basis functions for the representation of the incoming wavefront Φ and the measurements s = [s x , s y ] as derived in (8), the involved operators can be precomputed offline. These facts make the Landweber iteration coupled with a Kaczmarz strategy interesting for wavefront reconstruction from pyramid sensor data. A general convergence analysis of the linear Landweber-Kaczmarz method can be found in [45].
In the linear setting the Landweber-Kaczmarz method for wavefront reconstruction from pyramid wavefront sensor measurements (Algorithm 4) reads as:

Steepest descent-Kaczmarz approach
The idea of modified steepest descent algorithms coupled with a Kaczmarz strategy is comparable to the method described in [9] for non-linear problems. As in the previous method, we cyclically consider each measurement equation (5) and (6).
the steepest descent-Kaczmarz (SD-K) method for wavefront reconstruction using pyramid sensors (Algorithm 5a) is described by: Algorithm 5a Steepest descent-Kaczmarz method for pyramid sensors During an observation, the reconstructions have to be repeated up to 0.3 milliseconds. Assuming that the incoming wavefront do not change much from one time steps to the next and, in particular, tip & tilt do not change significantly, another idea (implemented in Algorithm 5b) would be to reconstruct in xdirection for even time steps t and in y-direction for the proximate odd time t+1 steps. The big advantage of Algorithm 5b consists in the reduction of the computational demand by more than 50% compared to the normal SD approach.

Algorithm 5b Modified steepest descent-Kaczmarz method for pyramid sensors if (t mod 2 = 0) do apply Algorithm 2 in x-direction only else if apply Algorithm 2 in y-direction only endif
The post-loop step of Algorithm 1-3, i.e., the averaging of the two reconstructions is not necessary for Algorithm 4, 5a, and 5b since we only obtain one reconstruction Φ K . Please note that for the Kaczmarz-type methods it is merely necessary to choose one initial guess Φ 0 instead of two as required for Algorithm 1-3.

Numerical results
We test the quality of the reconstruction approaches by continuously correcting the incoming wavefront in closed loop AO. In this setting, the wavefront sensor measures the incoming phase after passing the deformable mirror, i.e., the sensor sees the difference between the incoming wavefront and the correction induced by the mirror. For numerical simulations, we use the end-to-end simulation tool Octopus developed by ESO [48,49]. As already mentioned in the introduction, we test the performance of the proposed methods for an ELT-sized telescope system. In particular, we consider the METIS instrument on the 39 m sized ELT for non-modulated and modulated pyramid wavefront sensors having a 74 × 74 spatial sampling. Although the observing facility has a primary mirror diameter of 39 m, for METIS only the inner 37 m are used. The incoming wavefronts are simulated by a realization of the von Karman atmospheric model having 35 layers. The system runs at a frequency of 1 kHz for the non-modulated sensor and at a frequency of 500 Hz for the modulated sensor. The mirror geometry in the simulations corresponds to the M4 geometry planned for the ELT which was just recently incorporated in Octopus. For the temporal control of the algorithms we use a simple integrator and optimize the gains with a resolution of 0.1. As a quality measure we use the long-exposure (LE) Strehl ratio, which is computed as the average on-axis Strehl ratio for all performed time steps. The Strehl ratio is defined as the ratio of the peak aberrated image intensity from a point source compared to the maximum attainable intensity using an ideal optical system limited only by diffraction over the system's aperture. The maximum achievable value is 1.  The quality results of the algorithms are expressed in terms of long-exposure Strehl ratios at an observing wavelength of 2.2 µm (K-band). Note that according to the specifications of the METIS instrument, Kband is not included in the science range. Instead, observations are performed in L-band (at λ 1 = 3.0 µm, λ 2 = 3.7 µm), in M-band (at λ = 4.7 µm) and in N-band (at λ = 10.0 µm). For analysis purposes, however, we find it useful to have the output at a shorter wavelength as well. As such, we use λ = 2.2 µm in the K-band where the imaging is performed. In our numerical tests we evaluate the reconstruction quality in a range of photon flux levels between 50 and 10000 photons per subaperture per frame for median atmospheric conditions. The simulation parameters are summarized in Table 1. In order to speed up convergence to the closed loop, in the first 13 time steps we apply the CuReD reconstructor [65,66], which corrects mainly for the low frequencies in the wavefront.

Optimal step size choice for SD iteration in the context of WF reconstruction from pyramid data
Before we compare the reconstruction quality of all proposed methods, we investigate the optimal step size choice for the steepest descent algorithm applied to WF reconstruction. For that analysis we consider the METIS instrument on the ELT having a pyramid sensor without modulation incorporated. The simulation parameters are identical to those listed in Table 1. As photon flux, we use 10000 photons per subaperture per frame. The reconstruction quality is evaluated after 500 time steps using 5 SD-iterations for each reconstruction in order to find the optimal choice of the step size. As listed in Table 2, best results are obtained for the SD iteration combined with the classical steepest descent step size. The reason for the small number of performed iterations is (amongst others) related to the roof sensor approximation for modeling a pyramid sensor and discussed below in more detail.  Table 2: SD-reconstruction (Algorithm 2) results for the non-modulated sensor in the K-band after 500 time steps using different step sizes.

Simulated closed loop performance
Let us analyze the closed loop performance of the developed algorithms and compare their reconstruction quality. Our reconstruction methods are all based on a simplification of the full pyramid sensor model. As a consequence, after some iteration steps, the reconstructions suffer from an approximation error and depart from the true solution of the full pyramid sensor model although the residuals with respect to the simplified model continue to scale down during the iterations. Due to the fact that the full non-linear pyramid sensor model P consists of two terms P = P lin + P rest , where the first term again contains two terms P lin = R lin + S lin (see part I of the paper [36] for more details), a reduction of the roof sensor residual (15) can imply an error increase of S lin Φ + P rest Φ in the residual corresponding to the full pyramid sensor model Besides the approximation error, another error source, the data error, is present in the reconstruction process. It is inevitable to search for an adequate stopping criterion taking into account both the difference between the real pyramid sensor operator P providing the measurements s and the approximate operator R lin , which builds the foundation of the model-based reconstruction algorithms, as well as data errors. For choosing the regularization parameter in the generally non-linear problem of wavefront reconstruction from pyramid data, we discuss the usage of Morozov's discrepancy principle. Assume that the pyramid sensor provides noisy data s δ fulfilling s − s δ < δ for some noise level δ > 0. The iteration is terminated with stopping index k * (δ, s δ ) when for the first time the residual is below τ δ for some τ > 1, i.e., The discrepancy principle combined with a criterion for controlling the approximation error can be transferred to the application of only a few CGNE-or SD-iterations resulting in a very low value for k * as confirmed by a huge number of numerical simulations performed within this study. In particular, one iteration suffices to provide high reconstruction quality when using a warm restart of the system. That is, in the first time step, the initial guess is chosen as zero, i.e., Φ 0,0 := 0, and at time step t > 0 the initial value is set to the reconstructed phase of the previous step, i.e., Φ t,0 := Φ rec t−1 . By employing the reconstruction of the previous step as initial guess Φ 0 , we significantly decrease the computational complexity since Algorithm 1, 2 and 5 are scaled down to non-iterative gradient based methods by applying only one corresponding iteration step. The warm restart technique improves the convergence speed of the iterative solvers and additionally slightly increases the quality performance. The suitable number of iterations is also depending on the number of incident photons, since a high photon flux results in reduced data noise and vice versa.
In our applications a total number of K = 1 iterations turned out to be optimal with respect to the reconstruction quality and the computational complexity of the method. Except for the Landweber type approaches (Algorithm 3 and 4), we use more than one iteration, but already K = 5 Landweber steps  Table 3: Long-exposure Strehl ratios in the K-band obtained with the presented algorithms after 500 closed loop simulation steps for a pyramid sensor without modulation. Best results are obtained for the CGNE approach (Algorithm 1), the SD (Algorithm 2), and Landweber-Kaczmarz iteration (Algorithm 4).
combined with an adapted choice of the relaxation parameter and the warm restart technique are enough to obtain satisfying reconstruction quality.
In case of one CGNE-or SD-iteration the two algorithms coincide when using the classical steepest descent step size (12). Additionally, the step sizes in the SD method discussed in the previous Section do not differ for one SD-iterate except for the classical SD step size and the MG step size. Hence, for the numerical simulations with results provided in Table 3 and Table 4 we used the minimal gradient step size (13) in order to have an additional comparison of step size choices as well. Since Algorithm 5b has a reduced computational complexity compared to Algorithm 5a, we only consider the modified SD-Kaczmarz algorithm in our numerical tests. As above, one SD-Kaczmarz iteration suffices as well.
Numerical tests suggest that for METIS an interpolation to a finer grid than given by the subaperture spacing results in an increased reconstruction quality. In the XAO case the corresponding improvement was less significant. This may be related to the difference in subaperture sizes of both systems (21 cm in XAO versus 50 cm in METIS).
Corresponding results having a cold start (Φ 0 = 0) for every time step t can be found in Table 2 for Algorithm 2 utilizing a pyramid sensor without modulation while results using the warm restart technique for all presented algorithms are summarized in Table 3 for the non-modulated pyramid sensor and in Table 4 for the modulated pyramid sensor. Hence, the warm restart technique improves the reconstruction quality of the SD approach from an LE Strehl ratio of 0.8322 having a cold start to 0.8412 with the warm restart.
In case of zero modulation, best reconstruction quality is obtained for the Landweber-Kaczmarz approach using the two measurement sets alternating and for the gradient based approaches (Algorithms 1-2) calculating two reconstructions and averaging at the end. For the sensor having modulation 4 λ/D, surprisingly, the CGNE approach even outperforms the Landweber-Kaczmarz algorithm except for the simulations with 50 photons per subapertures per frame. However, the differences in the results are very small anyhow. In addition to the K-band results shown in Table 3, we provide the long-exposure Strehl ratios in other science bands as defined by the instrument specifications. Table 5 shows the quality in L-, M-, and N-bands obtained with the Landweber-Kaczmarz algorithm in the high flux case (10000 ph/subaperture/frame) for the non-modulated sensor.
The simulations for a modulated sensor whose results are presented in Table 4 were performed with a frame rate of 500 Hz. In order to have a direct comparison of the non-modulated and modulated sensor, we additionally run a simulation at a frame rate of 1 kHz (instead of 500 Hz) using a pyramid sensor with modulation 4 λ/D. For the modulated sensor with the CGNE method we obtain the LE Strehl ratio of 0.8782 in the K-band in the high flux case after 500 time steps and the LE Strehl ratio of 0.8415 for the non-modulated sensor. This result fits well our previous experiences with other model-based reconstruction algorithms according to which the modulated sensor provides a higher quality compared to the non-modulated one.
All in all, the developed reconstruction algorithms deliver comparable quality and allow for robust and accurate wavefront reconstruction with low computational costs.   Table 5: Long-exposure Strehl ratios in L-, M-, and N-bands obtained for the non-modulated pyramid sensor with the Landweber-Kaczmarz algorithm in the high flux case (10000 ph/subaperture/frame) after 500 closed loop simulation steps.

Comparison to interaction-matrix-based approaches
In the literature, there are many variants of interaction-matrix-based approaches: statistical estimators or least-squares methods; zonal or modal [47] control approaches (i.e., the degrees of freedom are modes or actuators/subapertures). The least-squares approach applied to the pyramid wavefront sensor reaches high correction accuracy without regularization, at least if the number of degrees of freedom is small as demonstrated in [25,24,23]. However, for large-scale AO systems the least-squares attempt turned out to be less accurate and the minimum variance estimator allowing statistical regularization is preferred [1,17] (see also MAP [11,28,47,50], MMSE [1]). One should note that regularization typically requires optimization (fine tuning) of the regularization parameters. Moreover, the regularized control matrix has to be recomputed each time the seeing conditions or the photon flux change, which is a rather time-consuming task. (As already mentioned, the computational complexity required for setting up the command matrix scales as O(N 3 ) and the application of this command matrix on the sensor data as O(N 2 )). Often, in practice and also in simulations, the non-modulated sensor being operated with an interactionmatrix-based approach, is reported to be unstable, see, e.g, [28,50]. One can for instance apply some tricks, like using a 'wrong' command matrix derived for the modulated sensor, or heavily fine-tune the regularization parameters to filter out the unstable modes in the correct interaction matrix (measured or computed for the sensor with modulation 0), which has to be performed on the fly and is a very time-consuming task.
Recently, there was a result published in [50] for the non-modulated sensor running in Octopus with a modal MVM at 1 kHz frame rate. The achieved quality in the K-band was reported to be 0.62 for the high flux case (10000 photons/subaperture/frame). For comparison, the pyramid sensor with modulation 4 λ/D was reported to provide in the same environment the LE Strehl ratio of 0.80. As recently reported in [32], another variant of MVM, the zonal minimum variance reconstructor in the YAO simulation tool [64], which is a zonal regularized approach, achieves LE Strehl of 0.89 in case of the modulated pyramid sensor (with modulation 4) and high photon flux.
Comparing the performances of the described algorithms, we can draw the following conclusions. For the pyramid sensor without modulation our reconstruction algorithms, which use the forward model of the sensor, allow not only to close the loop easily, but also to achieve a stable correction over time with a quality significantly higher compared to the interaction-matrix-based reconstructor utilized in [50]. In case of the pyramid sensor with the optimal amount of modulation, our algorithms achieve a reconstruction operation # of flops

Computational complexity
We define the computational complexity of the algorithm as a number of required floating point operations (flops). Let n denote the number of subapertures in one direction, then N = n 2 indicates approximately the number of unknowns to be found.

Complexity of Landweber iteration and Landweber-Kaczmarz iteration for pyramid sensors
We only consider the complexity of the operations that have to be performed online and exclude the pre-calculations needed in the application of the operators Q and Q * from our considerations. The number of floating point operations for every step in the Landweber iteration approach for wavefront reconstruction using pyramid sensors (Algorithm 3 and Algorithm 4) is provided in Table 6. The post loop step of the Landweber algorithm consists of finding the average between the two resulting reconstructions, which requires one summation and one division by a scalar. Altogether, this step is summed up to 2n 2 operations.
Since we perform the mentioned operations twice (in x-and in y-direction), for K iterations we obtain the complexity C 4 (n; K) = 8n 3 + 2n 2 · K for the Landweber-Kaczmarz approach (Algorithm 4) and C 3 (n; K) = 8n 3 + 2n 2 · K + 2n 2 flops for the application of the Landweber iteration (Algorithm 3) having the additional step of averaging.

Complexity of SD and SD-Kaczmarz algorithm for pyramid sensors
We again only consider the operations performed online. The complexity of one steepest descent iteration (Algorithm 2 and Algorithm 5) is indicated in Table 7. In case of the classical steepest descent iteration (Algorithm 2) a subsequent averaging Φ = 1/2(Φ x + Φ y ) has to be performed, which costs additionally 2n 2 flops. Therefore, the number of flops for the steepest descent-Kaczmarz approach applied to pyramid sensors (Algorithm 5a) is given by C 5a (n; K) = 12n 3 + 8n 2 + 2 · K, for the modified Algorithm 5b by C 5b (n; K) = 6n 3 + 4n 2 + 1 · K, and for the steepest descent approach (Algorithm 2) by C 2 (n; K) = 12n 3 + 8n 2 + 2 · K + 2n 2 , where K indicates the number of steepest descent steps.

Complexity of CGNE for pyramid sensors
The CGNE method consists of three steps: 1. a pre-computation and initialization step which have to be done for both x-and y-direction once, 2. the CG-loop for K iterations performed twice in x-and y-direction, 3. a post loop step in which we average the two obtained reconstructions.
The number of flops in the CGNE algorithm for pyramid sensors for operations which are not precomputed offline is indicated in Table 8. Summing up the specified operations for both s x and s y data, we see that the initialization step consists of 8n 3 flops, the loop of 8n 3 + 20n 2 + 4 · K and the post loop step of 2n 2 operations. Hence, the CGNE complexity for pyramid sensors sums up as C 1 (n; K) = 8n 3 + 20n 2 + 4 · K + 8n 3 + 2n 2 .
However, since the CG method is known to require the fewest number of iterations, K usually is smaller compared to, e.g., the Landweber iteration.
O N 3/2 C4(n; 5) 16, 3 · 10 6 325, 2 · 10 6 20 % not all actuators are actively controlled. However, for a theoretical comparison of complexities such assumptions are still relevant. Note that in principle C 1 (n; 1) = C 2 (n; 1) because the CGNE algorithm can already be terminated after the calculation of Φ 1 . However, we consider one full CGNE step in the Table. As a remark, we mention that in our comparison we have omitted the additional computational effort required in the presented model-based algorithms for computation of deformable mirror commands from the reconstructed wavefront shape. This step can be represented as a bilinear interpolation from the n × n grid of subapertures to the (n + 1) × (n + 1) grid of DM actuators, which requires 4(n + 1) 2 flops to be performed. Also, we would like to mention that the time-saving features of MVM approaches like parallelizability and pipelineability are valid in our algorithms as well.
The developed algorithms allow to significantly reduce the numerical effort of the wavefront reconstruction step in an AO loop compared to the computational load related to the solvers based on matrix-vector multiplication. This is illustrated especially well for the XAO system having a huge number of active actuators. The computational effort of MVM-based wavefront estimators is extremely demanding in this case. In contrast, the usage of analytically developed wavefront reconstructors allows one to heavily reduce the numerical effort of the AO loop. For instance as shown in Table 9, the modified steepest descent algorithm reduces the computational load of the wavefront reconstruction step in the XAO loop to approximately 3% of the MVM effort while still providing high reconstruction quality.

Conclusion and Outlook
In this paper we have studied the application of well-known iterative algorithms for solving the inverse problem of wavefront reconstruction from pyramid wavefront sensor data in the field of astronomical Adaptive Optics. From the performed end-to-end numerical simulations we can draw the conclusion that all studied algorithms deliver very similar reconstruction quality. However, it is preferable to apply the Kaczmarz versions of the algorithms or the CGNE approach, since they provide a slightly better reconstruction quality, though, the difference in the achieved quality between all the methods is minor. The best quality is obtained with the CGNE approach (Algorithm 1) and with the Landweber-Kaczmarz iteration (Algorithm 4), which at the same time is part of the slowest among the algorithms under comparison. If one decides to go for speed at the cost of a negligible quality loss, one should choose the modified steepest descent-Kaczmarz version combined with the classical step size choice (Algorithm 5b).
As shown by numerical results presented in this study, the proposed algorithms, which are partially iterative methods, allow to keep the numerical effort of the wavefront reconstruction step in an AO loop low compared to the computational load of solvers based on matrix-vector-multiplication. This has an especially big impact for the considered XAO system having a huge number of active actuators. For instance, the modified steepest descent algorithm reduces the computational load of the wavefront reconstruction step in the XAO loop to approximately 3% of the MVM effort while still providing high reconstruction quality.
Even when using simplifications of the pyramid sensor model, all proposed algorithms provide stable highquality reconstruction and (almost) reach the quality of interaction-matrix-based approaches in which the full pyramid model is assumed. Especially for the non-modulated sensor, the linear iterative algorithm give stable and very accurate wavefront reconstructions. If we compare the methods presented in this paper with the P-CuReD, all of them are outmatched by the P-CuReD with respect to both speed and quality. Nevertheless, in the proposed iterative methods, there is a possibility to investigate the full pyramid sensor model for future developments. Remarkable quality improvements are hoped for those adaptions. A big advantage of the iterative methods over the P-CuReD is that the full pyramid sensor model or real life features such as telescope spiders or the low wind effect can be incorporated. For the P-CuReD it may even be impossible to adapt the algorithm to a more sophisticated pyramid sensor model.
Finally, we would like to mention that investigations of the behavior of iterative algorithms in the presence of the so called "spiders" (support structures of the secondary mirror segmenting the telescope pupil into disjoint parts), the application of non-linear iterative algorithms for wavefront recnstruction as well as further quality evaluations to meet specifications of the METIS instrument are part of our further research [34,37,56].