## Abstract

We present two different implementations of the Fourier domain preconditioned conjugate gradient algorithm (FD-PCG) to efficiently solve the large structured linear systems that arise in optimal volume turbulence estimation, or tomography, for multi-conjugate adaptive optics (MCAO). We describe how to deal with several critical technical issues, including the cone coordinate transformation problem and sensor subaperture grid spacing. We also extend the FD-PCG approach to handle the deformable mirror fitting problem for MCAO.

©2006 Optical Society of America

## 1. Introduction

Multi-conjugate adaptive optics (MCAO) [1, 2] refers to the use of several deformable mirrors (DMs), each conjugate to a different altitude, to correct for volume atmospheric turbulence and thereby increase the corrected field of view in an astronomical adaptive optics system. By *wavefront reconstruction* we mean the determination of DM actuator commands, given wavefront sensor signals. In this paper we assume a linear model for sensor signals as a function of the turbulence profile, we assume DM displacements depend linearly on actuator commands, and we employ an optimal, or minimum variance, approach to wavefront reconstruction [3]. This gives rise to a pair of linear systems—one to estimate the volume turbulence profile and second to determine the actuator commands. Solution of the first system is referred to as the estimation step, or tomography, and solution of the second is called the fitting step.

With planning for a future 30-meter class telescope (TMT) that will rely on MCAO underway [4], it is important to be able to solve both the estimation step and the fitting step very quickly. Since TMT will require tens of thousands of sensors and actuators, conventional wavefront reconstructors, which are based on direct matrix-vector multiplication, may be prohibitively expensive. This has motivated the development of alternative approaches that rely on sparse direct matrix equation solvers [5] and on iterative methods like the preconditioned conjugate gradient algorithm (PCG) [6] with multigrid preconditioning (MG-PCG) [7, 8, 9] and Fourier domain preconditioning (FD-PCG) [10].

It should be noted that these iterative methods yield approximations to the optimal wavefront reconstruction. In simulations [9, 10] both MG-PCG and FD-PCG have been shown to be very rapidly convergent, so few iterations are needed. In addition, the cost per iteration is low, so both algorithms provide essentially optimal reconstructions at relatively low cost.

This paper is a follow-up to [10], where FD-PCG was first introduced. We deal with some important issues that either were inadequately addressed or not considered at all in the original paper. These include the fitting step, the cone-coordinate transformation problem, and problems related to computational grids that do not correspond to sensor subaperture corner locations. The cone-coordinate problem is especially important; failure to address it correctly will lead to field distortions in images obtained after MCAO correction.

This paper is organized as follows. In section 2 we review basic MCAO concepts and introduce notation. Section 3 contains a review of optimal volume turbulence estimation. In subsection 3.1 we present the cone-coordinate transformation problem. In subsections 4.1 and 4.2 we review the basic ideas underlying the FD-PCG algorithm and we present 2 different implementations of the algorithm. In subsection 4.3 we describe how to deal with computational grids that are finer than the sensor subaperture grid. Section 5 contains a review of the fitting step and an outline of how to efficiently apply FD-PCG to fitting. In subsection 5.1 we demonstrate FD-PCG fitting for a simulated MCAO system for a 30-meter class telescope. We close with a summary in section 6.

## 2. Basic Concepts and Notation

We assume a layered atmosphere with a refractive index profile which we denote by *ψ*(**x**,*z*_{l}
), *l* = 1,…,*n*
_{layer}, where the discrete layer heights *z*_{l}
are known, and **x** = (*x*,*y*) denotes lateral position. Since rapid refractive index variations are concentrated in regions of high turbulence, *ψ* is often referred to as the *atmospheric volume turbulence profile*. (Note that one can obtain the ${C}_{n}^{2}$ profile, a standard measure of turbulence strength as a function of height *z*, by integrating the square of *ψ*.) We further assume that at any instant in time, *ψ* is a realization of a wide-sense stationary stochastic process characterized by a Kolmogorov power spectral density function. See [11] for details.

The phase, or wavefront aberration, is obtained by propagating light from an idealized point light source, either at infinity in the case of a natural guidestar (NGS) or at a known finite altitude in the case of a laser guidestar (LGS), through the atmosphere to the pupil plane of the telescope. We denote dependence of phase on volume turbulence by

where *θ* denotes direction of the guidestar and **x** denotes location in the pupil plane. Using a geometric optics approximation, we model the propagator *P* as

In our model, both *ψ* and *ϕ* represent optical path differences and have units of length.

We assume that wavefront sensor measurements can be modelled as noisy, discrete approximations to the gradient of the phase. For example, suppressing the dependence on direction *θ*, the output of the (*i,j*)th element of a Shack-Hartmann wavefront sensor can be modelled [12] as

$${s}_{i,j}^{y}=\frac{1}{h}\left[\frac{(\varphi \left({x}_{i},{y}_{j}+h\right)-\varphi ({x}_{i},{y}_{j}))}{2}+\frac{\left(\varphi \left({x}_{i}+h,{y}_{j}+h\right)-\varphi ({x}_{i}+h,{y}_{j})\right)}{2}\right]+{\eta}_{i,j}^{y},$$

where *h* represents the sensor subaperture width and ${\eta}_{i\mathit{,}j}^{x}$
and ${\eta}_{i\mathit{,}j}^{y}$
represent sensor noise and grid discretization effects. Thus ${s}_{i,j}^{x}\approx \frac{\partial \varphi}{\partial x}({x}_{i},{y}_{j})$ and ${s}_{i,j}^{y}\approx \frac{\partial \varphi}{\partial y}({x}_{i},{y}_{j})$. We use the symbol Γ to represent the (discrete gradient) mapping from phase *ϕ* to sensor measurements **s** = (*s*^{x}
, *s*^{y}
).

Discretization in the direction variable *ψ* arises naturally because of the finite number of guidestars in a practical MCAO system. We impose a nodal discretization on *ψ* in the lateral variable **x** at each layer *z*_{l}
and nodal discretization of *ϕ* in each (discrete) direction *θ*. Hence the operators *P* and Γ have discrete matrix (approximate) representations. In addition, the Kolmogorov spectrum for (discretized) volume turbulence has a block diagonal matrix representation in the Fourier domain, with one diagonal block for each layer, due to independence of the layers.

In order to efficiently implement computational techniques based on the discrete Fourier transform, it is necessary to impose equispaced rectangular grids on the lateral variable **x** = (*x*,*y*) at each altitude *z*_{l}
. However, telescope apertures are typically circular or annular. Hence it is necessary to introduce an additional masking operator *M*. We take [*Ms*](*x*_{j}
,*y*_{j}
) to equal *s*(*x*_{i}
,*y*_{i}
) if (*x*_{i}
,*y*_{j}
) lies within the aperture and we take it to equal zero for (*x*_{i}
,*y*_{i}
) outside the telescope aperture. In an abuse of notation, we will also refer to the masking function whose value *M*(*x*_{j}
,*y*_{j}
) = 1 for (*x*_{j}
,*y*_{j}
) inside the aperture and *M*(*x*_{j}
,*y*_{j}
) = 0 for (*x*_{j}
,*y*_{j}
) outside the aperture.

The composition of discretized propagation, discrete gradient (corresponding to sensing), and masking will be represented by the volume turbulence-to-sensor measurement operator

By *atmospheric turbulence tomography*, or volume turbulence estimation, we mean the task of estimating the (discretized) volume turbulence profile *ψ* from noisy sensor measurements

In the implementation of the fitting step, one first propagates the light to the pupil plane from guidestars at infinity in several specified sample directions to obtain an estimated phase, *ϕ*
_{est}. The sampling must be dense enough to well-represent the phase, but not so dense as to make the computations overly expensive. The estimation plus this first stage of the fitting step can together be viewed as a glorified interpolation scheme to extend the phase to directions that have not been sampled by the guidestars.

The next stage of fitting is to determine deformable mirror (DM) actuator commands in order to compensate for the estimated phase. To this end, let *m*(**x**,*z*_{k}
) denote the displacement of the DM at conjugate altitude *z*_{k}
. The correction for pupil-plane phase aberrations due to the DMs can be represented as

Here *P* is again a geometric optics propagation operator with a representation as in Eq. (1); the DM displacements are analogous to turbulence layers placed at the DM conjugate altitudes. We assume the DM displacements can be represented as

The vector **a**
_{k} represents the actuator commands for the kth DM, and *H*_{k}
is the actuator-to-DM influence operator. We represent the combination of Eqs. (4) and (5) as

where **a** = (**a**
_{1},…,**a**
_{nDM}) is the concatenation of the actuator commands. As in the estimation step, nodal discretization of the DM displacements gives rise to a discrete matrix to represent *H*. The second stage of fitting can then be formulated as follows: Find the actuator command vector **a** for which *H*
**a** best matches (in a sense that we will precisely define in section 5) the estimated phase *ϕ*
_{est}.

## 3. Turbulence Estimation for an LGS-NGS MCAO System

To obtain sufficiently bright light sources to perform accurate volume turbulence profile estimation, it is necessary to make use of LGS beacons. Unfortunately, these beacons provide limited information about certain low-order components of the profile like tip-tilt, so additional NGS beacons are needed [5]. In this section we address technical issues that arise from a mixture of LGS and NGS beacons.

As in the previous section, we let *ψ* represent a nodal discretization of the volume turbulence profile. We separate wavefront sensor measurements into high-order components *s*_{h}
(from the LGSs) and low-order components *s*_{t}
(from the NGSs),

Each of *G*_{h}
, *G*_{t}
has form similar to *G* in Eq. (3). From [3], the optimal volume turbulence estimate is then given by

where ${\parallel x\parallel}_{A}^{2}\stackrel{\mathrm{def}}{=}{x}^{T}\mathit{Ax},{C}_{h}$ is the covariance matrix for high order sensor noise,*C*_{t}
is the covariance matrix for low order sensor noise, and *C*_{ψ}
is the covariance matrix for volume turbulence. The first two terms in (8) correspond to the pair of entries in (7), while the third is a *prior*, or regularization term, whose role is to incorporate prior statistical information about the volume turbulence.

By setting the gradient of the bracketed quantity in Eq. (8) equal to zero, we obtain

Following [5, 10], we express the covariance matrix for high-order noise in the presence of tip-tilt uncertainty as

where *N*_{h}
describes the noise within the high-order sensors, and *T* is a low-rank block matrix with 2*n*
_{LGS} columns (*n*
_{LGS} denotes number of LGSs). The product ${C}_{h}^{-1}$
*s*_{h}
can be interpreted as a noise-weighted, tip-tilt-removed version of the high-order sensor signals sh. We then rewrite Eq. (29) as

where

The two matrix products on the right hand side of Eq. (14) have rank 2*n*
_{LGS} and *n*_{t}
, respectively, where *n*_{t}
is the total number of low-order sensor measurements. Hence we can decompose

where *U*
_{1} has 2*n*
_{LGS} columns and *U*
_{2} has *n*_{t}
columns. On the other hand, the matrix *A*_{h}
in Eq. (13) has full rank. The Sherman-Morrison formula [6] allows to write

$$\phantom{\rule{3.5em}{0ex}}={A}_{h}^{-1}+\underset{L}{\underbrace{\left[{W}_{1}\phantom{\rule{.2em}{0ex}}{W}_{2}\right]{\left[\begin{array}{cc}I-{W}_{1}^{T}{U}_{1}& -{W}_{1}^{T}{U}_{2}\\ {W}_{2}^{T}{U}_{1}& I+{W}_{2}^{T}{U}_{2}\end{array}\right]}^{-1}\left[\begin{array}{c}{W}_{1}^{T}\\ -{W}_{2}^{T}\end{array}\right]}}$$

where

Then setting *b* = *b*_{h}
+ *b*
_{lr}, we obtain from (12) and (17)

From (17) we see that the product *Lb* can be computed by (i) taking a few dot products to obtain ${W}_{1}^{T}$
*b* and -*W*
^{2}
^{T}
*b*; (ii) applying the inverse of a small (2*n*
_{LGS} + *n*_{t}
× 2*n*
_{LGS} + *n*_{t}
) matrix to the result of (i); and (iii) adding together scalar multiples of the columns of *W*
_{1}, *W*
_{2}, where the scalars come from (ii). The product ${A}_{h}^{-1}$
*b* is computed iteratively using FD-PCG.

#### 3.1. The Cone-Coordinate Transformation

With geometric optics propagation, the contribution to the pupil-plane phase at position **x** due to light that has passed through a turbulence layer at altitude *z*_{l}
from a guidestar at height *H* with orientation *θ* is proportional to

Here *c* = 1 - *z*_{l}
/*H* and **s** = *z*_{l}*θ* are the cone compression factor and the shift, respectively, and the mapping

is referred to as the *cone coordinate transformation*.

From Eq. (20) we see that the cone coordinate transformation allows us to replace the combination of a rescaling and a shift with a simple shift. This simplifies the representation of the high-order geometric optics propagator Ph and facilitates the FD-PCG method. Figyre 1 illustrates the discrete grids that conform to cone coordinates and to standard coordinates.

The difficulty that arises with a (nodal) cone coordinate representation of the volume turbulence profile is that the low-order propagator *P*_{t}
then has a more complicated representation. To overcome this problem, we interpolate back to standard coordinates before applying *P*_{t}
. In cone coordinates, Eq. (8) takes the form

The *G̃*_{h}
in the first term on the right hand side denotes the cone-coordinate representation of *G*_{h}
in Eqn. (7). The *I*_{c}
in the second term denotes the transformation from cone coordinates back to standard coordinates. In practice all operators have discrete matrix representations and *I*_{c}
is a 2-D interpolation matrix, which is very sparse. In cone coordinates, Eqs. (12)-(15) must be slightly modified, with *G̃*_{h}
taking the place of *G*_{h}
and *G*_{t}*I*_{c}
taking the place of *G*_{t}
. Given *$\tilde{\psi}$ _{est}*, one computes

to get an estimate for the volume turbulence profile in standard coordinates.

## 4. FD-PCG Implementation of the Estimation Step

The effectiveness of the FD-PCG algorithm is a consequence of the fact that the components of the matrix *A*_{h}
in Eq. (13) nearly all have nice Fourier domain representations. In particular, the discrete gradient operator Γ corresponding to the Fried model (2) has a diagonal Fourier representation due to the Fourier shift theorem. The Fourier representer for the inverse Kolmogorov turbulence covariance matrix ${C}_{\psi}^{-1}$ is block diagonal with diagonal blocks. In cone coordinates, the high-order propagator *P*_{h}
has a Fourier representer which has a block decomposition with diagonal blocks, again as a consequence of the shift theorem. The mask *M* unfortunately does not have a compact representation, but it can be well-approximated by a scalar multiple of a diagonal matrix, as we demonstrated in [10].

We decompose the propagator *P*_{h}
into blocks *P*_{kj}
, where index *k* corresponds to LGS direction and index *j* corresponds to turbulence layer. Then from Eqs. (13) and (3), Ah also has a block decomposition,

with

the *δ*_{ij}
is the Kronacker delta, the regularization term *B*_{j}
corresponds to the inverse covariance matrix for the jth layer of turbulence, and Γ_{x}, Γ_{y} are the x- and y-components of the discrete gradient Γ.

In cone coordinates, each component propagator *P*_{kj}
is a simple shift operator. Hence it has a Fourier representation

where *F* represents the 2-D discrete Fourier transform on an *n*_{x}
× *n*_{x}
grid and *P̂*_{kj}
is a simple filter (component-wise scalar multiplication) operator. We assume the components Γ_{x}, Γ_{y} of the gradient operator have analogous Fourier-domain filter representations $\widehat{\Gamma}$
_{x}, Γ̂_{y}. In addition, we assume the tip-tilt removed noise covariance is scalar, ${N}_{h}^{-1}$ = ${\sigma}_{h}^{-2}$
*I*. Then

where *B̂*_{j}
= *FB*_{j}*F*
^{-1}, ∗ denotes complex conjugate transpose, and

$$\stackrel{\mathrm{def}}{=}{\tilde{S}}_{k}.$$

As in [10], we have approximated the Fourier transformed pupil mask *M̂* by a scalar multiple of the identity to obtain (29) from (28). For simplicity of presentation we have taken the scalar to be 1 in Eq. (29). We will revisit the issue of masking in Section 4.3 below.

#### 4.1. Direct Implementation of FD-PCG

PCG iteration to solve *A*_{h}*x* = *b* requires storage of 4 vectors, each the size of the right-hand-side *b*; see [6] for details. The dominant costs at each iteration are typically from the operator multiplications *h* = *A*_{h}*d* and the preconditioner applications *z* = *C*
^{-1}
*r*.

PCG iteration can be applied directly to compute *x* = ${A}_{h}^{-1}$*b* in Eq. (19). In this case, motivated by Eq. (27) and approximation (29), we take the preconditioner to be the matrix *C* having block components

Let *C* denote the matrix comprised of blocks *Ĉ*_{ij}
. We demonstrated in [10] that there exists a permutation, or reordering of rows and columns, for which *Ĉ* is block diagonal and the diagonal block size is quite small. (Diagonal block size will be addressed in detail in Section 4.3 below.) Hence *Ĉ* can be inverted directly and efficiently stored. The preconditioning step *z* = *C*
^{-1}
*r* is then implemented as follows:

- 2-D Fourier transforms are applied to the blocks of
*r*(note that blocks correspond to turbulence layers; there are*n*_{layer}blocks), yielding a vector*f*.The entries of

*r̂*are permuted, yielding a vector*r̃*. *r̃*is multiplied by the block diagonal inverse of the reordered*Ĉ*, yielding*ẑ*.The inverse permutation is applied to

*ẑ*, yielding*z*.- Inverse Fourier transforms are applied to the blocks of
*ẑ*to obtain*z*.

In the direct implementation of PCG, computations *h* = *A*_{h}*d* are carried out as sparse matrix-vector multiplications using the block decomposition (24). The propagators *P*_{kj}
correspond to simple shift-and-adds and need not be stored, and the matrices *S*_{k}
are very sparse. We employ Ellerbroek’s biharmonic approximation to the inverse turbulence covariance [5], so the regularization matrices *B*_{j}
are also very sparse. Computations *z* = *C*
^{-1}
*r* are dominated by the layer-wise 2-D Fourier transforms in steps (i) and (v).

#### 4.2. Transformed Implementation of FD-PCG and Comparison with Direct Implementation

Alternatively, one can apply PCG to the transformed system *Âx̂* = *b*, where *Â* has blocks given in (27) and b has blocks *b̂*_{j}
= *F*[*b*]_{j}. In this case the preconditioning step *ẑ* = *C*
^{-1}
*r̂* is carried out as in the direct PCG approach, but the Fourier transforms in steps (i) and (v) are omitted.

To carry out the step *ĥ* = *Âd̂*, we see from (27) that Fourier domain propagations can be carried out as component-wise product, or filter operations, involving (complex-valued) scalar multiplications. Similarly, the regularization terms *B̂*_{i}
give rise to filtering operations. However, we see from Eq. (28) that since the transformed pupil mask *M̂* has a full matrix representation, a Fourier transform / inverse Fourier transform pair is required for multiplication by each *Ŝ*_{k}
.

If *n*
_{layer} = *n*
_{LGS}, then the number of Fourier transforms applied during each PCG iteration is the same for the direct implementation and the transformed implementation. The cost of propagation and regularization is very similar using both approaches. However, the transformed implementation requires *n*
_{layer} additional Fourier transforms to obtain *b̂* from *b* and *n*
_{layer} additional inverse Fourier transforms to obtain *x* from *x̂*. Moreover, the transformed approach requires 4 complex-valued PCG vectors, while the direct approach requires 4 real-valued vectors, which need half as much storage.

#### 4.3. Grid Masking for High Order Sensor Subapertures

In our model for high-order sensors we assume sensor subapertures are square and have vertices (corners) which lie on the computational grid. To facilitate very accurate turbulence profile estimation, we take a computational grid that contains additional points that are not subaperture vertices. In order for the preconditioner *C* to closely approximate the operator *A*_{h}
, we must incorporate information about this grid structure in the Fourier domain. To this end we extend the subaperture vertex grid to the entire computational domain, as shown in Fig. 2, and we define the subaperture mask

We then replace the *S̃*_{k}
in Eq. (29) with

and we construct the preconditioner *C* from Eq. (30). As in [10], we reorder the rows and columns of the Fourier representer *Ĉ* according to spatial frequency, so the permuted *Ĉ* is block diagonal with diagonal blocks that are *n*_{b}
× *n*_{b}
, where *n*_{b}
= *n*
_{layer}(Δ*s*/Δ*x*)^{2}, with Δ*s* equal to the sensor subaperture spacing and Δ*x* equal to the spacing of the computational grid. For example, for a 6-layer atmospheric model with the grid shown in Fig. 2, Δ*s*/Δ*x* = 2 and the diagonal blocks of the permuted *Ĉ* are 24 × 24.

## 5. FD-PCG for the Fitting Step

Given a volume turbulence profile *ψ*
_{est}, one propagates to the pupil plane to obtain an estimated phase, *ϕ*
_{est}. Slightly modifying the notation introduced at the end of section 2, we model the masked DM corrections to the phase as

where *M* is a pupil mask, a is the actuator command vector, *H*
**a** represents PDM figures at the conjugate altitudes, and *P*
_{DM} represents propagation from the nsample virtual guidestars, through “layers” that correspond to DM displacements at the conjugate altitudes, to the pupil plane. *P*
_{DM} has an *n*
_{sample} × *n*
_{DM} block decomposition, the actuator influence operator *H* is block diagonal with *n*
_{DM} blocks, and *a* can be decomposed into *n*
_{DM} blocks.

The optimal actuator command vector is given by

$$\phantom{\rule{1.4em}{0ex}}={A}_{\mathrm{fit}}^{-1}\left({H}^{T}{P}_{\mathrm{DM}}^{T}\mathit{MW}{\varphi}^{\mathrm{est}}\right).$$

where

The purpose of the matrix *W* is to allow selective weighting in certain sampling directions, and the role of *R* is to stabilize the actuator commands. One can also accommodate fast tip-tilt mirrors or other low-order corrective elements. These give rise to low-rank terms and can be handled as in the estimation step using the Sherman-Morrison formula.

The matrix *A*
_{fit} takes the form of *A*_{h}
in Eq. (13) provided we identify *MP*
_{DM}
*H*, *W*, and *R* in (34) with *G*_{h}
, ${N}_{h}^{-1}$, and ${C}_{\psi}^{-1}$ in (13). In order to efficiently implement FD-PCG, components of *A*
_{fit} must be well-approximated by matrices with sparse Fourier domain representations. For this reason, we select *W* and *R* to have blocks with translation invariant structure (the matrix representers are then block circulant with circulant blocks, or BCCB). The propagator *P*
_{DM} should automatically be translation invariant, and many DM models assume that the blocks of *H* are translation invariant. As was the case in the estimation step, only the pupil mask does not have a sparse Fourier representer, but it can be approximated as before. We can again employ either of the two FD-PCG solutions strategies described in Sections 4.1 and 4.2.

## 5.1. An Illustrative Example

We present results for the fitting step applied to the simulated MCAO system for a 30-meter telescope described in [10]. This incorporates a 6-layer Cerro Pachon model for the atmosphere, 5 laser guidestars at 90 Km, 4 natural guidestars to deal with tip-tilt uncertainty, two DMs (one conjugate to ground and the second to 12 Km) with a corrected field of view of 30 arcseconds, and 1/2 meter sensor subaperture and actuator spacings. For the estimation step we now use 256 × 256 computational grids at each of the 6 layers, with twice the grid resolution in [10] (1/4 meter grid spacing on a square domain with a 64 meter width). The performance of the FD-PCG algorithm for estimation is much the same as in [10], so we present no additional details here.

We propagated the estimated volume turbulence profile to the pupil plane and then carried out the rest of the fitting step using the transformed implementation of the FD-PCG algorithm; see Section 4.2. The actuator influence functions, which determine the matrix *H* in Eq. (4), are taken to be bilinear splines with control points that coincide with the high-order sensor subaperture vertices (see Fig. 2). The matrix *H* has translation invariant structure and is quite sparse. FD-PCG performance is summarized in Fig. 3.

Let **a**
_{k} denote the approximation to **a**
_{opt} in Eq. (33) after k FD-PCG iterations, and let ${\varphi}_{k}^{\text{DM}}$ denote the corresponding DM correction to the phase, obtained by substituting **a**
_{k} for a in Eq. (32). By the *residual phase error* we mean *ϕ*
^{true} - ${\varphi}_{k}^{\text{DM}}$, where ftrue is the true phase, obtained by propagating through the (simulated) true volume turbulence profile to the pupil plane.

Note that at most 3 iterations are required to reduce the residual phase error to its asymptotic level, where the difference between **a**
_{k} and **a**
_{opt} is negligible and produces negligible change in ${\varphi}_{k}^{\text{DM}}$. Any further reduction in phase error would require adding more DMs or more DM actuators to reduce the fitting error or decreasing sensor noise to reduce the estimation error.

In [10] we observed similar behavior of the RMS phase error as we varied the number of FD-PCG iterations in the estimation step. The main difference was that 6 to 8 iterations were required to achieve asymptotic error levels for the estimation step, while only 2 or 3 iterations are required to reach asymptotic error levels in the fitting step.

## 5.2. Cost of Estimation and Fitting for the Example

For the example above, elements are required to store the nodal discretization of *ψ* in the estimation step. Standard implementations of PCG require 4 copies of *ψ*. There is some additional overhead to store the components of *A*_{h}
, but by taking advantage of special structure (sparsity and translation invariance), this can be kept relatively small. However, the *W*
_{1} and *W*
_{2} matrices in the low-rank matrix *L* in Eq. (17) require a great deal more storage. Each column of each *W*_{i}
has *N*
_{est} elements, and there are 2*n*
_{LGS} + *n*_{t}
= 2×5 + 8= 18 columns in *W*
_{1}, *W*
_{2} for the simulation. Hence the storage costs in the estimation step are dominated by the low-rank terms.

Again in the estimation step, the low-rank computation *L*_{b}
(which is O(*N*
_{est})) is much cheaper than the cost of computing ${A}_{h}^{-1}$
*b* via FD-PCG iteration. The FD-PCG costs are the cost per iteration multiplied by the total number of iterations. Dominant costs per iteration are multiplications by *A*_{h}
and by the inverse preconditioner *C*
^{-1}. Whether we use the direct approach in section 4.1 or the transformed approach in section 4.2 we must apply forward and inverse discrete Fourier transforms to block vectors with *N*
_{est} entries (each block is 256 × 256 in our simulation). Using the fast Fourier transform, the asymptotic cost for this is then O(*N*
_{est}log*N*
_{est}). There are some additional order *N*
_{est} costs, e.g., from the multiplication by the relatively small blocks of *Ĉ*
^{-1}.

While a detailed discussion of parallel implementation is beyond the scope of this paper, it should be noted that the layered structure of *ψ* and the block structure of the components of *A*_{h}
and of *Ĉ*
^{-1} provide ample opportunities to dramatically reduce computational costs in a parallel computing environment.

Storage requirements for the fitting step are significantly lower than for the estimation step. In the above simulation, there are only 2 DMs, so with a 256 × 256 nodal discretization of each DM, only ${N}_{\mathrm{fit}}\stackrel{\mathrm{def}}{=}6\times {256}^{2}\approx \mathrm{130,000}$ storage elements are needed to represent the DMs. PCG again requires 4 copies, but there are no low-rank terms.

Computational costs in the fitting step are also significantly smaller than for estimation. The asymptotic costs of the FFTs at eachfitting iteration are much smaller because *N*
_{fit} = *N*
_{est}/3 (due to 6 layers in the estimation vs 2 DMs in the fitting). In addition, the 2 or 3 FD-PCG iterations required for the fitting step are much fewer than the 6 to 8 iteration needed for estimation.

## 6. Discussion and Conclusions

In [10] we introduced the FD-PCG algorithm to efficiently carry out the estimation, or tomography, step in minimum variance wavefront reconstruction. Advances in this paper included

- A demonstration that the FD-PCG algorithm can be adapted to efficiently solve the fitting step in optimal MCAO wavefront reconstruction. We show that for a simulated MCAO system for a 30-meter telescope, at most 3 FD-PCG iterations are needed for the fitting step.
- Solution to the cone coordinate transformation problem. This problem arises with mixed LGS-NGS systems and can lead to field distortion if it is not handled correctly. Our solution involves a cone-coordinate-to-standard coordinate interpolation and is efficiently implemented with sparse matrix techniques.
- Solution of problems related to computational grids that have higher resolution than the high-order sensor subaperture grids. These problems are solved with a subaperture mask. This leads to a somewhat more complex preconditioner, but the additional computational overhead is small relative to the cost of the 2-D Fourier transforms.
- Two separate implementations of the FD-PCG algorithm. In the direct approach, PCG vectors lie in the spatial domain and are real-valued, while with the transformed approach the PCG vectors lie in the Fourier domain and are complex valued. The direct approach requires slightly fewer Fourier transforms and slightly less storage than does the transformed approach. Perhaps the only advantage of the transformed approach is that allows easy implementation of wave optics propagators. These have a sparse Fourier domain representation but have a convolution representation in the spatial domain.

## Acknowledgments

This research was supported in part by the Air Force Office of Scientific Research through grant F49620-02-1-0297 and by a grant from the Optical Sciences Company.

## References and links

**1. **J. M. Beckers, “Increasing the size of the isoplanatic patch with multi-conjugate adaptive optics,” in Proceedings of European Southern Observatory Conference and Workshop on Very Large Telescopes and Their Instrumentation, M. H. Ulrich, ed., Vol. 30 of ESO Conference and Workshop Proceedings (European Southern Observatory, Garching, Germany, 1988), pp. 693–703.

**2. **D. C. Johnston and B. M. Welsh, “Analysis of multi-conjugate adaptive optics,” J. Opt. Soc. Am. A **11**, 394–408 (1994). [CrossRef]

**3. **T. Fusco, J. M. Conan, G. Rousset, L. M. Mugnier, and V. Michau, “Optimal wave-front reconstruction strategies for multi-conjugate adaptive optics,” J. Opt. Soc. Am. A **18**, 2527–2538 (2001). [CrossRef]

**4. **R. G. Dekany, M. C. Britton, D. T. Gavel, B. L. Ellerbroek, G. Herriot, C. E. Max, and J-P. Veran, “Adaptive optics requirements definition for TMT,” Advancements in Adaptive Optics, edited by
D. B. Calia, B. L. Ellerbroek, and R. Ragazzoni, Proc. SPIE **5490**, 879–890 (2004). [CrossRef]

**5. **B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A , **19**, 1803–1816 (2002). [CrossRef]

**6. **G. Golub and C. VanLoan, *Matrix Computations, 2nd Edition*, Johns Hopkins University Press, 1989.

**7. **L. Gilles, C. R. Vogel, and B. L. Ellerbroek, “Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction,” J. Opt. Soc. Am. A , **19**, 1817–1822 (2002). [CrossRef]

**8. **L. Gilles, B. L. Ellerbroek, and C. R. Vogel, “Preconditioned conjugate gradient wave-front reconstructors for multi-conjugate adaptive optics,” Appl. Opt. **42**, 5233–5250 (2003). [CrossRef] [PubMed]

**9. **B. L. Ellerbroek, L. Gilles, and C. R. Vogel, “Numerical simulations of multi-conjugate adaptive optics wavefront reconstruction on giant telescopes,” Appl. Opt. **42** (2003), pp. 4811–4818. [CrossRef] [PubMed]

**10. **Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. **45**, No. 21 (2006). [CrossRef] [PubMed]

**11. **J. W. Hardy, *Adaptive Optics for Astronomical Telescopes*, Oxford University Press, 1998.

**12. **D. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. , Vol. **67** No. 3, (1977), pp.370–375 [CrossRef]