Local ensemble transform Kalman filter, a fast non-stationary control law for adaptive optics on ELTs: theoretical aspects and first simulation results

We propose a new algorithm for an adaptive optics system control law, based on the Linear Quadratic Gaussian approach and a Kalman Filter adaptation with localizations. It allows to handle non-stationary behaviors, to obtain performance close to the optimality defined with the residual phase variance minimization criterion, and to reduce the computational burden with an intrinsically parallel implementation on the Extremely Large Telescopes (ELTs).


Introduction
The prospect of ELTs (GMT [1],TMT [2], E-ELT [3]) with their related high-dimensional Adaptive Optics (AO) systems, such as eXtreme AO (XAO), MultiConjugate AO (MCAO) or even Multi-Object AO (MOAO), has aroused many developments in computationally efficient control algorithms.A good summary of the various techniques recently developed and their references can be found in [4].Linear Quadratic Gaussian (LQG) regulator is one of them.This control approach, though introduced a long time ago in AO [5], has truly found a resonance only during the last decade.This success is motivated by the optimal control solution provided by LQG, as well as its ability to account naturally for various perturbations such as vibrations [6][7][8], limitations such as saturations [9] or Deformable Mirror (DM) dynamics [10].
LQG has thus demonstrated improved performance compared to other control solutions, both in numerical simulations and experimental implementations.However, this control solution suffers from various limitations, the major one being the computational complexity especially in the non-stationary framework.Standard and brute force application of LQG is known to exhibit deterring computational cost, both in terms of off-line computations of the Kalman gain and on-line computations through multiple Matrix Vector Multiplications (MVM).Basically, assuming a number n of state variables in the system, complexity of the off-line computations is proportional to n 3 , and complexity of the on-line computations is proportional to n 2 .Moreover, this number n increases with the square of the telescope diameter in classical AO, and also with the number of reconstructed layers in tomographic approaches.As a consequence, reduction of computational complexity of Kalman Filter (KF)-based control solution has triggered much work in the recent years.The computational issue is significantly increased once considering non-stationary systems.As turbulence and system characteristics evolve with time, one would wish to update the control solution accordingly.For the LQG, this means redoing the off-line computations regularly, hence an additional computational cost.
To circumvent these limitations, the Ensemble Transform Kalman Filter (ETKF) has been previously introduced [11].This approach derives from methods developed in geophysics and meteorology in order to adapt the KF to large scale systems with non-stationary models.The ETKF does indeed provide a theoretical alternative to the KF for non-stationary systems but has also dramatic computational limitations for these large scale systems.Thus, a new control solution is discussed in which ETKF has been modified to implement a zonal and localized approach briefly proposed in [12].This present paper differs from the previous one in that we are giving a detailed description of this new approach and the required mathematical formalism with new simulations in the case of a 16 m telescope and the resolution of the differential pistons issue.The zonal approach is motivated by the sparse matrices involved as well as by the difficulties encountered when dealing with huge number of modes (typically, a few thousand to tens of thousands for ELT AO systems) such as Zernike or Karhunen-Loève modes (edge effects, numerical issue in modes computing).In addition, zonal approach also allows to use Fourier domain decomposition and WF reconstruction, benefiting from fast computation, though Fourier domain based approaches may suffer from edge effects on circular (annular) apertures, requiring particular handling [13,14].Localization here stands for decomposition of the system pupil into small domains on which local estimations are performed.All the local estimations are only combined in a final step, in order to deduce the full-pupil information.The result is thus a spatially distributed estimation based on the ETKF, leading to a hierarchical con-trol scheme.This approach allows a highly scalable implementation, as local computations are done with a reduced complexity and can be performed by a dedicated parallelized computation resource (CPUs/GPUs cluster).
This article is organized as follows.Section 2 recalls the main ideas behind using the LQG control and a KF for a classical AO system, its limitations and some solutions developed to overcome these drawbacks.Section 3 explains why the ETKF-based control law can handle non-stationary behaviors and can offer the optimality of the KF, but becomes unsuitable for computational reasons.In section 4 we present a new algorithm, the Local ETKF, and we explain why it is highly suitable for classical AO systems on ELTs (fast and naturally parallel implementation, non-stationarity).In order to assess the theoretical analysis and to demonstrate the potentiality of this new control law, some numerical simulations are detailed in section 5. Finally, conclusion and outlook of our work in progress are given in section 6.

The classic LQG control solution: an optimal control law for AO
It is now well-known that LQG provides the optimal control solution in AO, at least as long as the various underlying models (turbulence, AO system components, noise) are consistent with reality.To some extent, models can be approximated and LQG proves to be more efficient than any other control solution and to be also robust.Usually, it is still pointed out that this solution suffers from higher complexity and computational cost.Various works have been carried out to reduce this complexity and to propose faster computation even for large scale AO systems.As the Local ETKF basically derives from KF and uses similar formalism, in the following we first recall the basics of LQG and point out its main limitations.For the sake of simplicity, we discuss hereafter the LQG implementation in the framework of a classical AO system, although extension to tomographic systems has already been addressed [4,15].

Optimality criterion, discrete-time equivalence
Let us consider the simple AO system described in Fig. 1.The phases φ tur (t), φ cor (t) and φ res (t) represent respectively the incoming turbulent phase, the AO correction provided by the DM and the residual phase.In the astronomical framework considered here, the performance criterion consists in minimizing the variance of the residual phase (identified usually as the empirical variance computed over a sufficiently large amount of time), with respect to the DM controls u.Though turbulence is a continuous time phenomenon, the residual phase is usually integrated by the WaveFront Sensor (WFS) over a time period ∆T (frame), the AO control system is digital and correction is applied to the DM through a Zero Order Hold (ZOH) so that controls are constant over time period ∆T : see chronogram of the AO system in Fig. 2. In this context, Kulcsár has shown in [16] that the optimal turbulence correction by a sampled AO system is a problem that can be fully addressed in discrete-time.Defining for any continuous-time variable x(t) its discrete-time counterpart x k as: the optimal performance criterion comes down to minimizing the discrete-time criterion: This equivalence is independent from the control solution, and while this result has been demonstrated assuming infinite dynamics of the DM (instantaneous response) in [16], Correia has also extended it to the case of non negligible DM dynamics in [10].
As a conclusion, in the following we will consider the optimal criterion defined by Eq. ( 2), and we will consider for that a discrete-time control problem.We will thus assume that the WFS is linear and provides a measurement y k of the incoming wavefront φ k , integrated over the frame period ∆T , plus additional noise through the relation: where D is a linear operator accounting for the WFS measurement of the phase and w k is a discrete zero-mean Gaussian measurement white noise.For the sake of simplicity, we will also assume that the DM is linear, has no dynamics and can be described by its influence matrix N (the abscence of DM dynamics is clearly an optimistic hypothesis, though valid in most current AO systems: the impact of DM dynamics in the proposed control scheme is beyond the scope of this paper and should be addressed in the future).Its correction over a frame period ∆T is constant (due to the ZOH) and is therefore equal to: A two-frame delay AO system is thus considered as described in Fig. 2. Now, the AO control problem of minimizing the discrete-time criterion (Eq.( 2)) in order to provide the optimal AO correction of the incoming turbulent wavefront φ tur (t) finds its solution using the stochastic separation theorem [17].In a first step, the turbulent phase is estimated.The solution to this stochastic minimum-variance estimation problem finds its expression in a KF.Then, the estimated phase is projected onto the mirror's space with a simple least-square solution.The overall control scheme is thus a typical LQG regulator.

Mathematical formulation for a classical AO system
By using the structure of the LQG control and the associated KF described in [16,18], the turbulent phase ϕ tur can be defined on a zonal or a modal basis, and x k , the state vector at time k, contains 2 occurences of this turbulent phase and 2 occurences of the voltages u k : The stationary stochastic linear state-space model can be defined with these two equations: The first equation characterizes the dynamics of the turbulent phase described for instance by a first-order Auto-Regressive (AR 1) model (a higher order AR model could be also considered): The vector u k contains the voltages that are applied on the DM.The vector v k is a zero-mean white Gaussian model noise.The second equation is the observation equation: The vector y k contains the measurements of the residual phase ϕ res k = ϕ tur k − ϕ cor k .As we wrote in section 2.1, the estimation part can be separated from the control part.The optimal control u k can be obtained by separately solving a deterministic control problem and a stochastic minimum variance estimation problem.In Eq. ( 5), the state vector x k includes control voltages, but these components do not need estimation.Thus, in the following, the vector x k contains only the 2 occurences of the tubulent phase.The symbol k / k denotes analysis or update and the symbol k / k-1 denotes forecast or prediction.In the linear Gaussian case, the optimal solution of the update estimate of the turbulent phase ϕ tur is given by a KF: where ŷk/k−1 = C × xk/k−1 , and with the update estimation error covariance matrix: and the Kalman gain: ) Σ w is the measurement noise covariance matrix and C 1 is an extracted matrix from C, defined by C 1 = 0 D .The prediction estimate, at time k, is simply obtained with xk+1/k = A 1 xk/k where A 1 is an extracted matrix from A, defined by A 1 = A tur 0 Id 0 .The prediction estimation error covariance matrix is calulated through the Discrete Algebraic Riccati matrix Equation (DARE): where Σ v is the model noise covariance matrix.The deterministic control problem is simply solved through a least-squares projection of the predicted phase φtur k+1/k (upper part of xk+1/k ) onto the DM's space.With a basic KF implementation, we must solve the DARE in order to calculate this Kalman gain H k [16,18].As we defined a stationary model, all the matrices of the state-space model are time-constant during the observation and the asymptotic formulation of the KF is applied without loss of optimality: the Kalman gain is then precalculated with an offline computation by solving the DARE, which gives the asymptotic Kalman gain (H k = H ∞ ).

LQG limitations for AO systems on ELTs
LQG has proved to provide efficient control and improved performance compared to other control solutions in the numerical simulations as well as in the experimental validations.Since the very first works in AO [5], various developments have been done in classical AO [16,19,20] and in wide field tomographic AO [18,19,21,22].Experimental validations in laboratory have been also carried out both in AO [6,18] and MCAO [15,23], while first validations on sky [7,24,25] provided good demonstration of performance.
Nevetherless this control solution suffers from various limitations, the major one being the computational complexity, especially in the non-stationary framework.

Numerical cost and computational complexity with stationary models
As recalled in introduction, brute force implementation of LQG is prohibitive in terms of computational complexity.As a consequence, reduction of computational cost of KF-based control solution has drawn much work in the recent years.For instance, Poyneer [26] has proposed the use of Fourier decomposition of atmospheric turbulence and the statistical independence of the Fourier modes to define a mode-by-mode regulator.A modal KF is used.This approach benefits naturally from mode decoupling (sparsity) and possible parallelization leading to computation speed up.Correia [27] improved the off-line computation of the Kalman gain and real time operations by replacing MVM by spectral iterative methods with sparse approximation of the turbulence covariance matrix which avoids the resolution of the DARE.Massioni [28] proposed the Distributed Kalman Filter (DKF), which approximates the Kalman gain by assuming the telescope pupil as the cropped version of an infinite-sized phase screen.The method is based on a Fourier transform that is used to decompose an infinite-order system into an infinite set of low, finite-order ones.Fourier transform is used as an off-line tool to accelerate Kalman gain computation, but on-line computation is also improved by sparsity in a zonal representation of phase.Gilles compares in [4] the performance and costs of two computationally efficient Fourier based tomographic wavefront reconstruction algorithms, the DKF and the iterative Fourier Domain Preconditioned Conjugate Gradient combined with a Pseudo-Open-Loop Control, showing very good performance of DKF and an acceptable computational burden.At the end, most of these solutions tends to propose a drastic reduction of the numerical complexity which becomes O(n × log(n)).

Dealing with non-stationary turbulence
LQG (but also control approaches that do rely on turbulence statistics models) is usually or systematically derived in a time-invariant framework: the turbulence model and the system models (WFS, DM responses) are considered as time-invariant.The turbulent model is usually predefined and the model matrices are usually derived from information gathered on the spot right before observation [7], or deduced from global statistics.This leads obviously to significant simplifications, both in formalism and computation: in particular, one can compute an asymptotic Kalman gain, as mentioned before, which significantly reduces the real-time complexity of LQG.This is also motivated by a pragmatic trade-off between accuracy and efficiency: the goal is to provide good performance with simple enough models of turbulence, thus usually mentioned as control-oriented models.Solutions mentioned in the previous section in order to speed up off-line and on-line computations rely on a time-invariant framework.
While we can expect system models to present slow time evolution (or at least predictible), this can not be said about turbulence models as turbulence statistics (seeing conditions, average wind profile ...) clearly evolve over time.Then, one can consider updating the turbulence models, once in a while, based either on dedicated instruments information (seeing monitor) or real-time data as it is already performed on the SPHERE AO system for Tip-tilt control [24,29].
In the latter case, closed-loop data are directly used to identify periodically the turbulent Tip-tilt models as well as possible vibrations.The resulting identification is used to update the Kalman models and control solution.However, this implies a full update of turbulence models and recomputation of the Kalman gain, as well as correct management of transitions [25].While this is affordable on limited-dimensional systems such as SPHERE, considering large scale AO systems comes to a dead-end.Recent developments of fast computation of Kalman gain [27,28] may provide an interesting way out.
However time-invariant system is not a prerequisite, and time-variant LQG can also be derived similarly: in the second case, the turbulent model matrices are time-variant, and consequently the Kalman gain computation can not be precomputed off-line by convergence of the DARE (Eq.( 12)).Therefore, it must be computed at each step based on the current values of the turbulent model and system matrices through Eqs. ( 11) and (12).Of course this situation leads to an unsuperable computational complexity as well as to the problem of defining, on every other steps, these time-evolving matrices.Thus, embedding turbulence models identification and update within the control scheme, without significant increase of computational cost while ensuring performance, still represents an unreachable Graal.
The Local ETKF-based solution aims both at proposing a time-variant KF-based control solution and a reduced complexity, benefiting from the intrinsically parallel algorithm.

The ETKF: large scale systems adaptation and non-stationary models possibility
Although the LQG approach for AO systems with a stationary model has been successfully implemented on 4-8 m class telescopes, a similar implementation with non-stationary models is not possible for complex AO systems on ELTs.The major reason is the transition to very high-dimensional systems when explicit storage and manipulation with theoretical covariance matrices are not possible.In the KF-based method, there are two difficulties in inverting 11).The first one is the size: it is a p × p matrix (p is the number of measurements).For complex AO systems on ELTs, it can be costly to compute the inverse when p = O (10 5 ).The second one is the fact that this matrix will be ill-conditioned and it may be extremely difficult to accurately evaluate its inverse.This situation led us to adapt a new method for large scale AO systems with non-stationary models.

Main ideas and mathematical formulation
The original Ensemble Kalman Filter (EnKF) method is based on the KF's Eqs. ( 9) and (10), except that the Kalman gain (Eq.( 11)) is calculated from the covariance matrices provided by an ensemble of model states: it is a Monte Carlo method [30,31].A fixed number m of model states of the turbulent phase on the whole pupil of the telescope composes the members of the ensemble X = {x 1 ; ...; x m }, and by integrating it foward in time, it is possible to calculate the empirical estimation error covariance matrices with the following statistical estimator: where x i is the mean of the ensemble's members and Z is called the matrix of the anomalies.To be precise, Z k/k−1 is the prediction anomalies matrix and Z k/k is the update anomalies matrix.In other words, the prediction estimate x k/k−1 (respectively the update estimate x k/k ) of the turbulent phase is given by the mean of the m members of the ensemble X k/k−1 (respectively X k/k ).All these members will constitute together a cloud of points in the state space and the spreading of this ensemble characterizes the prediction error variances by ).However, in this EnKF formulation, it has been shown the need to add in the measurements y k some perturbations (random variables calculated from a distribution with a covariance matrice equal to Σ w : for more details see [30]), but this use of 'pertubated measurements' introduces sampling errors that reduce the update covariance matrix accuracy.
Another variant has been developed in order to form the update ensemble X k/k deterministically, avoiding part of the sampling noise.In particular, the update estimation is not given by the mean x k/k of the m members (model states) of this ensemble X k/k , but directly by the estimate xk/k from the KF Eq. ( 9).In the following, we present this deterministic algorithm for transforming the prediction ensemble X k/k−1 into an update ensemble X k/k : it is called the Ensemble Transform Kalman Filter (ETKF).

The update step
The ETKF-based method doesn't explicitly calculate the empirical covariance matrices, but transforms prediction anomalies Z k/k−1 into update anomalies Z k/k by an Ensemble Transform Matrix (ETM) T k with the relation: such that the empirical prediction covariance matrix Σ k/k−1 and the empirical update covariance matrix Σ k/k , defined by Eq. ( 13), both match the theoretical Eqs. ( 10) and (11).There are different solutions for the ETM T k and following the mathematical approach in [32][33][34][35], a general form that satisfies Eqs. ( 10) and ( 11), and Eqs.( 13) and ( 14), is: We emphasize that, in this study, we assume the measurement noise covariance matrix Σ w to be a stricly positive diagonal matrix (no correlation between different subapertures).Given the Eigen Value Decomposition (EVD) of the matrix Id whose size is m × m, the solution for the ETM T k can be obtained with this relation: where the orthogonal matrix Q k contains the normalized eigenvectors of the m × m matrix in the square brackets of Eq. ( 15) and the diagonal matrix Γ k contains the eigenvalues of the EVD.Since Id ) is a positive definite real symmetric matrix, its eigenvalues are real, strictly positive and its eigenvectors are orthogonal.
Both prediction and update covariance matrices belong to the same linear subspace spanned by the ensemble prediction anomalies Z k/k−1 .A distinguishing feature of the update anomalies Z k/k produced by the ETKF is that they are orthogonal under the inner product (defining also a Euclidian norm): However, we have to notice that no more than m − 1 independent update anomalies are generated from the expression in Eq. ( 13) because the sum of the m prediction anomalies (columns of Z k/k−1 ) is equal to zero: therefore, the rank of the empirical covariance matrices Actually, in this update scheme with the ETKF-based control law, the Eq. ( 9) of the update estimate xk/k with the Kalman gain H k given by Eq. ( 11), is completely transformed by using the anomalies matrix Z k/k−1 and the ETM T k with the Sherman-Morrison-Woodbury identity (see Appendix A).As Σ w is a strictly positive diagonal matrix, it is straightfoward to calculate Σ −1/2 w .Let us define the vector and the matrix the expression for the update estimate is therefore: The matrix of the m members of the update ensemble X k/k is obtained with:

The prediction step
In order to obtain each of the m members of the prediction ensemble, we have to compute: where v i k+1 is a zero-mean random Gaussian vector characterising the model noise (with a covariance matrice equal to Σ v ).The prediction estimate x k+1/k is given by the mean of the m members of the prediction ensemble X k+1/k and the prediction anomalies matrix Z k+1/k is still given by the expression of Eq. ( 13).

Theoretical numerical complexity with non-stationary models
Let us note n, the dimension of the state vector x k and p, the dimension of the observation vector y k in Eq. ( 9).In order to determine the theoretical numerical cost with a non-stationary turbulence model, we have to calculate the total number of all the MVM computed during the update step and the prediction step.Actually, the significant computational cost of the ETKF method is the update step.It can be proved (Appendix B) that the resultant theoretical cost is: which is therefore linear over n and p with a proportional factor equal to m 2 .It will actually depend on the relative magnitudes of the parameters m, n and p.With this ETKF-based method, it is numerically more efficient if the value of m remains smaller than the values of n and p.

ETKF limitations for AO systems on ELTs
In the ETKF-based method, there are two related limitations [36] by using an ensemble with m members for the calculations of empirical covariance matrices with Eqs. ( 13) and (14).
The first limitation is the model space dimension of the finite size ensemble X much lower than the one on the whole pupil of the telescope.If the ensemble X has m members, then the empirical prediction covariance matrices Σ k/k−1 describe uncertainty only in the (m − 1)dimensional subspace spanned by the ensemble.The global update will allow adjustments to the system state only in this subspace which is usually rather small compared to the total dimension of the model space on the whole pupil of the telescope.As a small ensemble has few degrees of freedom available to represent estimation errors, this situation leads to sampling errors and a loss of accuracy with an underestimation of the true covariance matrices.Thus, these empirical covariance matrices calculated from the m members will not be able to match the model rank (the effective number of degrees of freedom of the model) unless this number of members increases considerably: this leads to an extremely high theoretical numerical complexity and no possibility of a realistic implementation on ELTs.The idea is therefore to perform locally all the updates so that, with different linear combinations of the ensemble members in various domains, the global update explores a much higher dimensional space.
The second (though related) limitation is the spurious correlations over long spatial distances produced by a limited size ensemble X .The empirical covariance matrices are calculated with the statistical estimator (Eq.( 13)) where the ensemble X of model states is considered as a statistical ensemble and the models (turbulence, WFS ...) are imperfect.Therefore, the correlation calculations from the ensemble sample assign non-zero values to correlations between variables separated by a large distance (compared to the Fried parameter r 0 ), which leads again to an underestimation of the true covariance matrices.It can be shown that variances of these spurious random correlations decrease when the number of members increases: however, it is not suitable for an implementation on ELTs.The idea is to determine the system's characteristic correlation distance, and then perform again locally so that the local update should ignore ensemble correlations for distances larger than this correlation distance/length (which is not equal to r 0 ).

The Local ETKF : an intrinsically parallel algorithm
In the previous section, we have described two related shortcomings with the ETKF-based control law, which lead to underestimate the theoretical covariance matrices.In order to overcome these drawbacks, we have proposed the necessity to use local updates.Thus, a new version called the Local ETKF [34,37] has been developed using domain decomposition and localizations.There are two common localization methods in the EnKF-based approaches.The first one is the Local Analysis (LA): the local update is performed explicitly by considering only the measurements from a local region surrounding the local domain by building a virtual local spatial window.This is equivalent to setting ensemble anomalies outside a local window to zero during the update.The second one is the Covariance Localization (CL) where the local update is performed implicitly: the prediction covariance matrices are tapered by a Schur-product with a distance-based correlation matrix.It increases the rank of the modified covariance matrix and masks spurious correlations between distant state vector elements.As the two methods yield very similar results [38], we present in this section the LA: thus, the idea of the Local ETKF is to split up the pupil of the telescope into various local domains on which all calculations of the update step are performed independently.

Description of the update step
In order to understand the principles of the Local ETKF, let us take the following example with a 16 m diameter telescope: Fig. 3 shows the locations of the valid actuators (blue dots) on the pupil sampled by a 32×32 SH-WFS with a Fried geometry: therefore, each square area between 4 valid actuators represents a subaperture.In the following, we describe the Local ETKF-based control law on a zonal basis: the turbulent phase is estimated on the locations of all DM's actuators.In this example, the domain decomposition is made of 25 (red) estimation subpartitions, there is an independent update.Therefore, this update step scheme enables highly parallel computations of markedly less data.We emphasize that it is essential to avoid as far as possible discontinuities between two neighboring local updates coming from two different local domains.To be precise, the two results of two update estimates coming from two nearby actuators on the border of two different local domains have to be similar.This can be ensured by choosing similar sets of observations for neighboring actuators, i.e. by taking two overlapping observation sets coming from two different neighboring local regions.Moreover it is necessary to have a smoothed localization by gradually increasing the uncertainty assigned to the observations until beyond a certain distance they have infinite uncertainty and therefore no influence.This can be done by multiplying (in Eq. ( 22) hereafter) the square root inverse of the measurement noise covariance matrix by a decreasing function (from one to zero) as the distance of the observations from the center of the local observation region increases [37,38].

Basic mathematical formulation of the Local ETKF
We introduce a new notation where a tilde denotes the local version of the corresponding variable.For example: Z k/k−1 means the local prediction anomalies used for a given local domain d; y k means the local measurements from the local observation region including and surrounding a given local domain d.In order to facilitate the understanding, we give also the size of each local variable obtained in the calculations.Therefore, we use the new variables n d and p d , which are defined as n and p: they are respectively the dimension of the local state vector and the dimension of the corresponding local measurement vector on a local domain d.
Thus, for the update step, on each local domain d, we have to determine independently each local update xk/k by computing: and With the concatenation of all those small matrices X k/k , we can finally compute the update members matrix X k/k .Then, the prediction step is computed globally on the whole pupil (as it is done with the ETKF in section 3.1.2):the computation of the m predicition members (Eq.( 19)) for the prediction matrix X k+1/k , the prediction estimate x k+1/k and the prediction anomalies matrix Z k+1/k .
For the sake of understanding the ability of this algorithm to handle non-stationary behaviors, each local update estimate (Eq.( 24)) can be rewritten by factoring S inov (to the right): where H k is called the local Kalman gain on the local domain d.With this last expression, we can clearly see that all local Kalman gains are computed during each update step, without the need to resolve a DARE.Therefore, we can have a non-stationary command with small modifications of the matrices Z k/k−1 and Σ w : there will be far fewer transitional issues than what happens with the LQG command and the full recomputation of all control matrices.Moreover, during the prediction step, by using a (to be defined) identification procedure, we can update the AR model (matrices A 1 and Σ v ) and take into account this modification with Eq. ( 19).

Theoretical numerical complexity with non-stationarity models
Let us note n max the largest value of the dimensions n d of all local update estimates on all local domains, and p max the largest value of the dimensions p d of all observation vectors on all observation regions.With the Local ETKF, we have to distinguish two computational costs.The theoretical cost of the update step on each local domain is: This complexity is again linear over n max and p max with a proportional factor equal to m 2 .We emphasize that the values of n max and p max can be significantly smaller than n and p when the pupil of the telescope has been split up into many subpartitions: therefore, by using the result given in Appendix B.1, the total number of operations on each local domain can be much smaller.Moreover, it will only depend on the size of the subpartitions (number of actuators) and not on the diameter of the telescope (see table 1 in section 5.3).Actually this total number of operations will depend on the relative magnitudes of m, n max and p max which are determined by an optimal trade-off between the size of the local domains, the size of the observation regions and the expected image quality (see the various simulations in section 5.2).
For the prediction step, the theoretical cost is: A key point is that both the calculations during the update step and the prediction step can be easely parallelized, which considerably speeds up the algorithm for AO systems on ELTs (see discussion in section 5.3).

Local phase estimation and differential pistons issue
One particular consequence of the Local ETKF approach is the differential pistons issue.Indeed, as described in section 4.1, we perform a partition of the pupil of the telescope into various local domains, on which one local update estimate of the turbulence phase will be produced.Afterwards, one global update estimate (for each member of the ensemble) and only one global prediction estimate are performed on the whole pupil.Turbulence continuity from a local domain to its neighbors, is ensured to some extent in the overlapping of the observation regions, that includes for each connected local domain, the surrounding subapertures.However this continuity does not extend to piston which is not measured by the WFS.
In other words, turbulence estimation on local domains is piston-free, which means ipso facto that a piston continuity issue will arise when combining all local update estimates of the turbulent phase.In order to calculate these differential pistons and to remove them, we have now implemented a first effective and fast algorithm inside the AO loop, using a least-squares based method, which globally minimizes the local discontinuities and ensure phase continuity on the whole pupil.
Though straightforward, this algorithm allows improving significantly the Local ETKF performance (see section 5.1).However, this algorithm does not yet take into account the noise propagation at the level of turbulence phase estimation, and thus the influences on local pistons.Of course, this is a foreseen improvement, that could simply rely on a regularized least-square method, taking advantage of our knowledge of phase estimation accuracy, embedded in the empirical update estimation error covariance matrix given by Σ k/k .

First simulation results with a classical AO system (D = 16 m)
In order to validate some of the items discussed in sections 3 and 4, we have implemented the three previous control laws (based on the KF, the ETKF and the Local ETKF) under the OOMAO Matlab environment [39] and we have also developed our OpenMP parallelized version.For this Single Conjugate AO (SCAO) system simulation, we consider a 16 m diameter telescope (without central obscuration) with a 32×32 microlenses S-H WFS.Only the subapertures with a surface enlightened more than 50 % are considered and the number of valid subapertures is 812 that means p = 1624.The linear response of the S-H WFS is emulated by a matrix which calculates differences of the phase in two directions at the edges of each subaperture [18,28].The AO system works in a closed-loop at 500 Hz.There is a two-frame delay between measurement and correction.We assume that the DM has an instantaneous response and the coupling factor of the actuators is 0.3.Using a zonal basis, the phase is estimated only on the actuators' locations and the number of valid actuators is 877 that means n = 1754.
The criterion used in the comparisons is the coherent energy, defined by E coh = exp(−σ 2 res ), where σ 2 res is the temporal mean of the spatial variances of the residual phase ϕ res .The loss of performance (in %) with the Local ETKF is calculated by using the optimal solution given by the KF: For the simulation of the atmosphere, we consider a Von Karman turbulence with a stationary model: r 0 = 0.525 m (at 1.654 µm), L 0 = 25 m, λ = 1.654 µm (for both WFS's and observation's wavelengths).Using Taylor's hypothesis, we can generate a superimposition of 3 turbulent phase screen layers moving at 7.5 ms −1 , 12.5 ms −1 and 15 ms −1 , with a relative C 2 n profile of 0.5, 0.17 and 0.33 respectively.Phase screens are generated on a 320×320 grid, with 10×10 points per each subaperture: the correction phase is therefore obtained by multiplying the prediction estimate on the actuators' locations with an influence matrix of the DM.
For the turbulence temporal model in the KF, in the ETKF and in the Local ETKF-based control laws, we have chosen a first order Auto-Regressive model (AR1) [16,18].Each value of the coherent energy has been calculated with several simulations of 5000 iterations (10 sec).

Convergences of the ETKF and the Local ETKF to the KF
The mathematical convergence of the Ensemble Kalman Filters to the Kalman Filter in the limit for large ensembles has been proved in [40] for linear forecast models.Therefore, the first step is to obtain this theoretical result with the simulations for the ETKF: it was done with one given value of noise variance equal to 0.04 Rad 2 for which the KF gives a coherent energy equal to 88.5 %.In Fig. 4, each curve gives the coherent energy loss in % between the performance obtained with a control law (ETKF or Local ETKF on a 25 domains partition without removing or by removing differential pistons) and the one given by the KF.This loss depends on the number of members and decreases to zero when the value of m increases.The performance provided by the Local ETKF without removing differential pistons (red dashed curve) is better than the one provided by the ETKF until a given number of members (≃ 626).Indeed, concerning the phase estimation with the Local ETKF, there is a problem of reattachment between the various local domains discussed in section 4.4: as the local update estimates of the turbulent phase are done separately on each local domain, we have to take into account the differential pistons and to reconstruct the entire phase estimate on the whole pupil of the telescope with all these various small turbulent phase estimates.
As the red solid curve shows, once differential piston has been removed with our first leastsquares based method, the coherent energy loss compared with the KF is always lower than 4 % for numbers of members greater than 100.With only 122 member, this loss is already about 3.5 % on a 25 domains partition (see also Fig. 5).

Performance of the Local ETKF with different partitions of the pupil of the telescope
The second step is to study different kinds of partitions and their influences on the performances for a fixed value of noise variance.Figure 5 shows the coherent energy losses (compared with the KF) obtained with the Local ETKF in the case of five kinds of partitions and various values of members in the ensemble.As already shown in Fig. 4, for a given partition, the larger is the number of members, the smaller is the loss of performance.In the same way, for a fixed value of members, the larger is the number of subpartitions (or the smaller is the number of actuators per local domain), the smaller is the loss of performance: it confirms the specific problems of the ETKF explained in section 3.3 and the solution given by the Local ETKF.The localization improves the efficiency of the ETKF-based approach because, firstly the ensemble needs only to encompass the uncertainty within each of the small local domains (with small numbers of degrees of freedom), and secondly the spurious long-range correlations produced by a limited ensemble size are removed.What is also important in this simulation with a 16 m telescope, is that we can choose a partition of the pupil, for which the loss of performance is less than 3%: for example with the 9 × 9 partition and a number of members equal to 101 (it will be also true with a larger number of subpartitions or a larger number of members).The third step is to study the influence on the performance for different values of noise variance with a fixed number of members (m = 101).Let us recall that in order to remove the differential pistons, our first least-squares based method doesn't take into account the measurement error covariance.Of course, for a given partition, when the noise variance increases, the loss of performance increases too.Nevertheless, as Fig. 6 shows, the larger is the number of subpartitions, the smaller is the loss of performance when the noise variance increases.Moreover, for the 9 × 9 partition (with a maximum of 16 actuators per local domain), the rosbustness of the Local ETKF compared to the KF seems to be very good for a large range of noise variance, while the correction of the differential pistons can be still improved.
Compared performance along the time for the KF and the Local ETKF is proposed in Fig. 7 for various kinds of partitions.This figure shows that convergence speed and stability of the Local ETKF is very similar to the KF.It underlines once more that increasing the number of subpartitions allows to reduce the loss of performance compared with the KF and that a very small difference in performance is brought going from a 9 × 9 to a 17 × 17 partitions.

First speed-tests for a parallel implementation
The complexity given in section 4.3 is only theoretical.We want to evaluate the performance of a real parallel implementation of this intrinsically parallel algorithm based on the Local ETKF: therefore, we made some runtime tests in order to demonstrate that this method could be used for a 40 m telescope with a frequency up to 500 Hz using nowadays computers.
Let us consider the 16 m telescope for which the loss of coherent energy compared with the KF is less than 3 % with the 9 × 9 partition (total of 69 domains) and 101 members.In the case of a 40 m telescope with a 80 × 80 SH-WFS (and the same r 0 ), by keeping local domains with the same maximum number of actuators, we can reach a similar coherent energy loss with the same number of members (m = 101).For the 16 m telescope with the 9 × 9 partition, there is no more than 16 actuators per domain.For the 40 m telescope, in order to have no more than 16 actuators per domain, we must take a 21 × 21 partition (total of 373 domains).In this AO configuration, our numerical simulations have confirmed that the loss of coherent energy is indeed less than 3%.Each cycle of the AO loop has 3 steps: all the local updates (calculated independently on each local domain), the prediction and the calculation of voltages u k .Therefore, we propose to consider a multi-core architecture (single computer or computation cluster), where each core is assigned to only one local domain and performs the local update computations.By using the results of Table 3 and Table 4 in Appendix B, Table 1 gives the theoretical numbers of floating-point operations during one cycle of the AO loop.Table 2 gives real speed-tests with a first version of our OpenMP parallelized code: we used a workstation with two Intel(R) Xeon(R) E5-2680 v2 CPUs (total of 20 cores).We have measured the time required for calculating, the update step for one local domain by using only one core, the prediction step by using all the 20 cores and the voltages u k by using also all the 20 cores.In the AO configuration on a 16 m telescope, we have only 69 domains.Nowadays we can have a single computer with 120 cores (8×15-cores CPUs) similar to the cores in our test workstation.Table 2 shows that, with this kind of single computer (which has more cores than the number of requested domains), our algorithm needs a total of 1.18 ms for one cycle of the AO loop.Therefore, we can easely implement and use our algorithm with a 500 Hz frequency for the AO loop.For the AO configuration on a 40 m telescope, the situation is not basically different: there are many more local domains, but during the update step, computations on each local domain take exactly the same time as in the 16 m case.The time required for the prediction step and the u k calculation became not negligible: however both of these computations (prediction estimate x k+1/k and voltages u k ) are rather well scalable with a larger number of cores.Right now, we cannot have a single computer with so many cores (373), and we have to consider computation cluster which implies additional time for communication stage between the update step and the prediction step.On our computation cluster with infiniBand 4x QDR network, this communication stage takes ∼ 2 msec.However the last generation of infiniBand (12x EDR) is claimed to be almost 10 times faster, which makes times for communication reasonable small and a real implementation on a 40 m telescope achievable.

Conclusion
We have presented a new control law, the Local ETKF, based on the LQG approach, a KF adaptation with localizations for large-scale AO systems and the assumption of a perfect DM dynamics.The advantages of this proposed method are significant.Firstly, as all the local Kalman gains can be calculated with empirical covariance matrices at each cycle of the AO loop, it enables to deal with non-stationary behaviors (turbulence, vibrations).Secondly, as the structure of this algorithm is intrinsically parallel, its implementation on ELTs can easily be done on a parallel architecture (CPUs/GPUs cluster) with a reduced computational cost: let us remind that the complexity of the update step does not depend on the diameter of the telescope but only on the maximum number of actuators per each local domain, which significantly speeds up the algorithm.The more subpartitions in the pupil of the telescope, the better is the performance, both in terms of coherent energy and of runtime.In our simulations with a Von Karman turbulence and a SCAO system with an adequate partition of the pupil, the loss of coherent energy compared with the optimal solution given by the KF can be less than 3 %.We have presented numerical simulations in the case of a SCAO system with an AR1 turbulence model on a zonal basis, but we can already extend this method with an AR2 model.For the runtime tests we have already developed an OpenMP parallelized version, but we are also currently working on a new version with GPUs on the COMPASS platform [41].
In the short-term, there are different points to study more in details.The first one is the influence of the turbulence charateristics on the size of the local observation regions which reflects the distance, called the localization length, over which the correlations calculated with the ensemble's members are not meaningful.In the same way, we need to evaluate the influence of the spider arms on the choice of the subpartitions on the pupil of the telescope.In order to resolve the problem of differential pistons, we have to improve our fast least-squares based method by taking into account the error covariance matrices.Then, the robustness of the Local ETKF must be also studied when the parameters of the turbulence change.We must indeed estimate the impact of the turbulence model error (wind speed, turbulence profile, L 0 ), first in absence of turbulence model update.Then, assuming some turbulence model identification and update, one shall consider the gain brought by this non-stationary control solution, its stability and the speed of convergence.
In the long-term, we have of course to demonstrate the potentials of the Local ETKF in the case of wide field tomographic AO and to consider the extension to limited DM dynamics.Afterwards, two other important aspects already developed in geophysics must be also explored in AO.The first one is when the operator C in Eq. ( 6) is non-linear [42].The EnKF-based method does not require a linearized model, an advantage over the KF-based method, and this could be very suitable for non-linear WFS.The second one is the possibility of asynchronous observations assimilation [37,43], for the multi-rate case in the prospect of Natural Guide Star and Laser Guide Star wavefront sensing for wide field tomographic AO.
Actually the significant computational expense is the EVD and the 3 matrix-matrix multiplications: Total number of multiplications: (m 2 + m) × n + (m 2 + 7m + 41) × p + 2m 3 + 3m 2 + (3 + γ)m The number of additions is the same order of magnitude as the number of multiplications.

B.2 The Prediction step
By using a zonal basis with an AR1 turbulence model, A 1 is composed by two n act × n act diagonal blocks, one of them is the identity matrix.Therefore, multiplying A 1 with a vector is reduced to only one multiplication with the block A tur (the other one consists on a copy).The voltage u k is calculated with the MVM: u k = (N T N) −1 N T × φtur k+1/k .The matrix P = (N T N) −1 N T is a sparse matrix on a zonal basis.In our SCAO system, for a 16 m (respectively a 40 m) telescope, the sparsity of this matrix is 51 % (respectively 88 %).Actually, on a zonal basis, the matrix P can be much more sparse when, for each actuator, we take into account only the neighboring actuators close to less than 2 pitches.

Fig. 3 .
Fig. 3. Partition of the actuators domains on a 16 m telescope pupil with a 32×32 SH-WFS.

Fig. 4 .
Fig. 4. Convergences of the ETKF and the Local ETKF performances to the KF.

Fig. 5 .
Fig. 5. Coherent energy loss with different partitions and a fixed noise variance.

Table 1 .
Numbers of operations with m = 101 during one cycle of the AO loop.

Table 2 .
Runtime (in msec) with m = 101 during one cycle of the AO loop.

Table 3 .
Numbers of multiplications during the update step.