Distributed Kalman filtering compared to Fourier domain
preconditioned conjugate gradient
for laser guide star tomography on
extremely large
telescopes
Luc Gilles, Paolo Massioni, Caroline Kulcsár, Henri-François Raynaud, and Brent Ellerbroek
Luc Gilles, Paolo Massioni, Caroline Kulcsár, Henri-François Raynaud, and Brent Ellerbroek, "Distributed Kalman filtering compared to Fourier domain preconditioned conjugate gradient for laser guide star tomography on extremely large telescopes," J. Opt. Soc. Am. A 30, 898-909 (2013)
This paper discusses the performance and cost of two computationally efficient
Fourier-based tomographic wavefront
reconstruction algorithms for wide-field laser guide star (LGS) adaptive optics
(AO). The first algorithm is the iterative Fourier domain preconditioned
conjugate gradient (FDPCG) algorithm developed by Yang et al.
[Appl. Opt. 45, 5281
(2006)], combined with pseudo-open-loop
control (POLC). FDPCG’s computational cost is proportional to N
log(N), where N denotes the dimensionality
of the tomography problem. The second algorithm is the distributed Kalman filter
(DKF) developed by Massioni et al. [J. Opt. Soc. Am. A 28, 2298
(2011)], which is a noniterative spatially
invariant controller. When implemented in the Fourier domain, DKF’s cost is also
proportional to N log(N). Both algorithms are
capable of estimating spatial frequency components of the residual phase beyond
the wavefront sensor (WFS) cutoff frequency thanks to regularization, thereby
reducing WFS spatial aliasing at the expense of more computations. We present
performance and cost analyses for the LGS multiconjugate AO system under design
for the Thirty Meter Telescope, as well as DKF’s sensitivity to uncertainties in
wind profile prior information. We found that, provided the wind profile is
known to better than 10% wind speed accuracy and 20 deg wind direction accuracy,
DKF, despite its spatial invariance assumptions, delivers a significantly
reduced wavefront error compared to the static FDPCG minimum variance estimator
combined with POLC. Due to its nonsequential nature and high degree of
parallelism, DKF is particularly well suited for real-time
implementation on inexpensive
off-the-shelf graphics processing units.
You do not have subscription access to this journal. Cited by links are available to subscribers only. You may subscribe either as an Optica member, or as an authorized user of your institution.
You do not have subscription access to this journal. Figure files are available to subscribers only. You may subscribe either as an Optica member, or as an authorized user of your institution.
You do not have subscription access to this journal. Article tables are available to subscribers only. You may subscribe either as an Optica member, or as an authorized user of your institution.
You do not have subscription access to this journal. Equations are available to subscribers only. You may subscribe either as an Optica member, or as an authorized user of your institution.
Goal is to compute turb. est
and subsequently
from
5
global tip/tilt removal with gradient
weights
6
tomography right-hand-side vector
7
initial guess is previous frame
estimate
8
initial residual
(-mult)
interpolation offsets
in unit
of mesh size and measured from
lower left nearest
neighboring grid pt
,
has irregular
stencil (18 weights per subap)
Laplacian has regular
stencil (five nonzero weights)
13
preconditioning step
14
()
2D forward FFTs for each layer
15
read
values from
read 4 pts on each os layer and 1 pt on
each non-os layer
(see Fig. 1)
16
apply
MVM
17
write
values
()
back to
write result back to original
locations
18
()
2D backward FFTs for each layer
19
initial search direction
Table 2.
FDPCG Tomography Algorithm Iterations and Output Approximate Solution
20
For
loop of FDPCG iterations
21
imaged search direction
(-mult)
22
dot product
23
dot product
24
solution update
25
If
if test on iteration number
26
residual update
27
preconditioning step
28
dot product
29
search direction update
30
End
End
end of if test on iteration number; end
of for loop on iteration
number
31
solution after
iterations
32
optional two-step prediction using wind
profile estimate (bilinear interpolation)
Table 3.
DM Fitting and Temporal Filtering
33
ray tracing along
fitting
directions
(bilinear interpolation)
34
fitting right-hand-side vector
35
unfiltered DM command
36
filtered DM command (LPF)
Table 4.
DKF Tomography Algorithm Description
88
DKF Computations
Goal is to compute turb. est
and subsequently
from
and from
89
measurement signal from tomography
estimate
90
pseudo-open-loop error signal
91
global tip/tilt removal with gradient
weights
92
Kalman gain mult
93
SD convolutions on 2D arrays
94
convolutions via 2D FFTs and 2D FD
filter component-wise
array mults
95
solution update
96
one-step prediction using wind profile
estimate to obtain current frame
turbulence
estimate
to be stored in memory
97
one-step prediction using wind profile
est and ray tracing
along
directions
98
[B33–B36]
fitting right-hand-side vector,
unfiltered command and final
temporally
filtered command
Table 5.
Required Atmospheric and Aperture-Plane Grid Sizes for LGS MCAO
Observations with TMT at 60 deg Zenith Angle and for a 2 arcmin
Diameter Fitting Field of Viewa
37
B
C
D
E
38
1D nb of phase
points per layer
at 60 deg
zenith angle for 2 arcmin diameter fitting field (square
boxes enclosing metapupils)
FD3 or DKFP
FD3o2
FD3o6 or DKFo6P
39
(layer 1)
61
123
123
40
(layer 2)
71
143
143
41
(layer 3)
80
80
161
42
(layer 4)
90
90
179
43
(layer 5)
110
110
221
44
(layer 6)
122
122
245
45
(aperture-plane grid for )
123
46
Mem 1e6 real numbers
Such grid sizes are large enough to accommodate LGSs as low as
80 km in altitude.
Table 6.
Computation Requirements for Three FDPCG Iterations for the Case
,
2, ,
Denoted, Respectively, as FD3, FD3o2, and FD3o6a
47
FDPCG Computations
48
FD3
FD3o2
FD3o6
FD3o6 Formulas
49
Nlgs
6
Nb. of LGSs
50
Nps
6
Nb. of estimated phase screens
51
Nos
0
2
6
Nb. of layers oversampled
52
Ndm
2
Nb. of DMs
53
Ng/1e6
0.03
Nb. of gradients
54
nnz(W)/Ng
2.00
Noise weighting
55
Cost W mult
0.06
56
20.18
DM interaction matrix (bicubic
IF)
57
Cost Ga mult
1.25
58
Mem Ga (Compressed Sparse Row)
9.66
59
Niter
3
Nb. of CG iterations
60
N/1e6
0.05
0.08
0.20
61
Nfft
128
mixed
256
1D nb. of phase pts on FFT
grids
62
Cost of single Nfft x Nfft FFT
0.23
NA
1.05
63
nnz(Gamma)/Ng
8.07
aperture-plane WFS grad matrix
64
Mem Gamma
0.95
65
N0/1e6
0.02
66
3.39
3.12
2.93
ray tracing LGS-to-aperture
67
2.10
1.95
1.84
68
Cost open-loop grads
[B3]
1.28
69
Cost tip/tilt bias removal
[B5]
0.06
70
Cost tomo rhs vec
[B6]
2.16
2.01
1.91
71
nnz(L)/N
4.96
4.97
4.97
Laplacian matrix
(regularization)
72
Cost L mult
0.25
0.38
0.38
73
Cost tomo matrix mult/iter
[B8]
4.75
4.73
4.52
74
Cost precond FFTs/iter [B14,
B18]
2.75
6.03
12.58
75
Cost precond MVM/iter
[B16]
2.36
9.44
37.75
76
Cost precond/iter
[B13]
5.11
15.47
50.33
77
Cost tomo matrix mult
[B8,B21]
19.00
18.92
18.06
78
Cost precond
[B13,B27]
15.34
46.40
150.99
79
Mem precond
2.63
9.75
37.50
80
Cost dot prods, vec adds
[B8,B22–B24,B26,B28,B29]
0.80
1.23
1.21
81
Cost total tomo
38.64
69.91
173.52
82
Nfit
9
Nb. of fiting directions
83
3.70
3.49
3.49
ray tracing fitting
field-to-aper
84
Cost H mult [B33]
3.02
2.85
2.85
85
Prefitting Cost grand total
41.66
72.76
176.38
86
Quadratic Incr. from FD3o6 (nm) with
no wind
profile
knowledge ()
40.00
20.00
0.00
87
Quadratic Incr. WFE from FD3o6 (nm)
with perfect
wind profile
knowledge
35.00
10.00
Negative sign indicates reduced
WFE
The number of operations is given in units of MMACs per frame,
and memory requirements in units of megabytes (MB). The last
column contains the formulas used for FD3o6 (letters refer to
columns and numbers to rows).
Table 7.
DKF Tomography Algorithm Computational Cost Analysis for the Case of
the Kalman Gain Multiplication Implemented via FFTs for the Case
,
,
Denoted, Respectively, as DKFP and DKFo6Pa
99
DKF Comps with Kalman Gain Mult Implemented via
FFTs
100
DKFP
DKFo6P
DKFo6P Formulas
101
N/1e6
0.05
0.20
E60
102
Nfft
128
256
103
N0/1e6
0.02
C65
104
Cost open-loop grads
[B3]
1.28
C68
105
Cost tomo grads
[B89]
2.10
1.84
E67
106
Cost grad est. error
[B90]
3.41
3.16
107
Cost tip/tilt bias removal
[B91]
0.06
C69
108
Nb. of convolution kernels
72
109
Cost of Kalman gain mult FFTs
[B94]
4.13
18.87
110
Cost of Kalman gain mult FD filter
mults [B94]
4.72
18.87
(complex-valued mults)
111
Total Cost Kalman gain mult
8.85
37.75
112
Mem FD Kalman gain
9.00
36.00
(complex-valued)
113
Cost phase vector add
[B95]
0.10
0.40
tomo estimate
114
Cost shift matrix mult
[B96,B97]
0.40
1.62
predictions
115
Cost total tomo
14.00
43.87
116
Cost H mult [B33]
3.02
2.85
E84
117
Prefitting Cost grand total
17.02
46.72
118
Quadratic Incr. WFE from FD3o6 (nm)
with perfect wind
profile
knowledge
18.00
Negative sign indictated reduced
WFE
119
Quadratic Incr. WFE from FD3o6 (nm)
with no wind profile
knowledge
()
60.00
42.00
120
Quadratic Incr. WFE from FD3o6 (nm)
with 20 deg wind
direction
error and 10% wind speed error
41.00
The last column contains the formulas used for DKFo6P (letters
refer to columns and numbers to rows).
Table 8.
DKF Tomography Algorithm Computational Cost Analysis for the Case of
the Kalman Gain
Multiplication
Implemented via SD Convolutions, for the Case
,
,
Denoted, Respectively,
as DKFP-SD and
DKFo6P-SDa
121
DKF Comps with Kalman Gain Mult Implemented via
SD Convolutions
122
DKFP-SD
DKFo6P-SD
DKFo6P-SD Formulas
123
nK
60
120
1D nb of points in convolution
kernels
124
Ng/1e6
0.03
0.12
125
Cost grad est. error
[B90]
3.41
3.16
D106
126
Cost tip/tilt bias removal
[B91]
0.06
C69
127
Cost SD Kalman gain mult
[B93]
669.60
10713.60
128
Cost SD Kalman gain mult/Cost FD
Kalman gain mult
75.68
283.81
D127/D111
129
Mem SD Kalman gain
0.99
3.96
The last column contains the formulas used for DKFo6P-SD (letters
refer to columns and numbers to rows).
Tables (8)
Table 1.
FDPCG Tomography Algorithm Initialization Steps
4
FDPCG Computations
Goal is to compute turb. est
and subsequently
from
5
global tip/tilt removal with gradient
weights
6
tomography right-hand-side vector
7
initial guess is previous frame
estimate
8
initial residual
(-mult)
interpolation offsets
in unit
of mesh size and measured from
lower left nearest
neighboring grid pt
,
has irregular
stencil (18 weights per subap)
Laplacian has regular
stencil (five nonzero weights)
13
preconditioning step
14
()
2D forward FFTs for each layer
15
read
values from
read 4 pts on each os layer and 1 pt on
each non-os layer
(see Fig. 1)
16
apply
MVM
17
write
values
()
back to
write result back to original
locations
18
()
2D backward FFTs for each layer
19
initial search direction
Table 2.
FDPCG Tomography Algorithm Iterations and Output Approximate Solution
20
For
loop of FDPCG iterations
21
imaged search direction
(-mult)
22
dot product
23
dot product
24
solution update
25
If
if test on iteration number
26
residual update
27
preconditioning step
28
dot product
29
search direction update
30
End
End
end of if test on iteration number; end
of for loop on iteration
number
31
solution after
iterations
32
optional two-step prediction using wind
profile estimate (bilinear interpolation)
Table 3.
DM Fitting and Temporal Filtering
33
ray tracing along
fitting
directions
(bilinear interpolation)
34
fitting right-hand-side vector
35
unfiltered DM command
36
filtered DM command (LPF)
Table 4.
DKF Tomography Algorithm Description
88
DKF Computations
Goal is to compute turb. est
and subsequently
from
and from
89
measurement signal from tomography
estimate
90
pseudo-open-loop error signal
91
global tip/tilt removal with gradient
weights
92
Kalman gain mult
93
SD convolutions on 2D arrays
94
convolutions via 2D FFTs and 2D FD
filter component-wise
array mults
95
solution update
96
one-step prediction using wind profile
estimate to obtain current frame
turbulence
estimate
to be stored in memory
97
one-step prediction using wind profile
est and ray tracing
along
directions
98
[B33–B36]
fitting right-hand-side vector,
unfiltered command and final
temporally
filtered command
Table 5.
Required Atmospheric and Aperture-Plane Grid Sizes for LGS MCAO
Observations with TMT at 60 deg Zenith Angle and for a 2 arcmin
Diameter Fitting Field of Viewa
37
B
C
D
E
38
1D nb of phase
points per layer
at 60 deg
zenith angle for 2 arcmin diameter fitting field (square
boxes enclosing metapupils)
FD3 or DKFP
FD3o2
FD3o6 or DKFo6P
39
(layer 1)
61
123
123
40
(layer 2)
71
143
143
41
(layer 3)
80
80
161
42
(layer 4)
90
90
179
43
(layer 5)
110
110
221
44
(layer 6)
122
122
245
45
(aperture-plane grid for )
123
46
Mem 1e6 real numbers
Such grid sizes are large enough to accommodate LGSs as low as
80 km in altitude.
Table 6.
Computation Requirements for Three FDPCG Iterations for the Case
,
2, ,
Denoted, Respectively, as FD3, FD3o2, and FD3o6a
47
FDPCG Computations
48
FD3
FD3o2
FD3o6
FD3o6 Formulas
49
Nlgs
6
Nb. of LGSs
50
Nps
6
Nb. of estimated phase screens
51
Nos
0
2
6
Nb. of layers oversampled
52
Ndm
2
Nb. of DMs
53
Ng/1e6
0.03
Nb. of gradients
54
nnz(W)/Ng
2.00
Noise weighting
55
Cost W mult
0.06
56
20.18
DM interaction matrix (bicubic
IF)
57
Cost Ga mult
1.25
58
Mem Ga (Compressed Sparse Row)
9.66
59
Niter
3
Nb. of CG iterations
60
N/1e6
0.05
0.08
0.20
61
Nfft
128
mixed
256
1D nb. of phase pts on FFT
grids
62
Cost of single Nfft x Nfft FFT
0.23
NA
1.05
63
nnz(Gamma)/Ng
8.07
aperture-plane WFS grad matrix
64
Mem Gamma
0.95
65
N0/1e6
0.02
66
3.39
3.12
2.93
ray tracing LGS-to-aperture
67
2.10
1.95
1.84
68
Cost open-loop grads
[B3]
1.28
69
Cost tip/tilt bias removal
[B5]
0.06
70
Cost tomo rhs vec
[B6]
2.16
2.01
1.91
71
nnz(L)/N
4.96
4.97
4.97
Laplacian matrix
(regularization)
72
Cost L mult
0.25
0.38
0.38
73
Cost tomo matrix mult/iter
[B8]
4.75
4.73
4.52
74
Cost precond FFTs/iter [B14,
B18]
2.75
6.03
12.58
75
Cost precond MVM/iter
[B16]
2.36
9.44
37.75
76
Cost precond/iter
[B13]
5.11
15.47
50.33
77
Cost tomo matrix mult
[B8,B21]
19.00
18.92
18.06
78
Cost precond
[B13,B27]
15.34
46.40
150.99
79
Mem precond
2.63
9.75
37.50
80
Cost dot prods, vec adds
[B8,B22–B24,B26,B28,B29]
0.80
1.23
1.21
81
Cost total tomo
38.64
69.91
173.52
82
Nfit
9
Nb. of fiting directions
83
3.70
3.49
3.49
ray tracing fitting
field-to-aper
84
Cost H mult [B33]
3.02
2.85
2.85
85
Prefitting Cost grand total
41.66
72.76
176.38
86
Quadratic Incr. from FD3o6 (nm) with
no wind
profile
knowledge ()
40.00
20.00
0.00
87
Quadratic Incr. WFE from FD3o6 (nm)
with perfect
wind profile
knowledge
35.00
10.00
Negative sign indicates reduced
WFE
The number of operations is given in units of MMACs per frame,
and memory requirements in units of megabytes (MB). The last
column contains the formulas used for FD3o6 (letters refer to
columns and numbers to rows).
Table 7.
DKF Tomography Algorithm Computational Cost Analysis for the Case of
the Kalman Gain Multiplication Implemented via FFTs for the Case
,
,
Denoted, Respectively, as DKFP and DKFo6Pa
99
DKF Comps with Kalman Gain Mult Implemented via
FFTs
100
DKFP
DKFo6P
DKFo6P Formulas
101
N/1e6
0.05
0.20
E60
102
Nfft
128
256
103
N0/1e6
0.02
C65
104
Cost open-loop grads
[B3]
1.28
C68
105
Cost tomo grads
[B89]
2.10
1.84
E67
106
Cost grad est. error
[B90]
3.41
3.16
107
Cost tip/tilt bias removal
[B91]
0.06
C69
108
Nb. of convolution kernels
72
109
Cost of Kalman gain mult FFTs
[B94]
4.13
18.87
110
Cost of Kalman gain mult FD filter
mults [B94]
4.72
18.87
(complex-valued mults)
111
Total Cost Kalman gain mult
8.85
37.75
112
Mem FD Kalman gain
9.00
36.00
(complex-valued)
113
Cost phase vector add
[B95]
0.10
0.40
tomo estimate
114
Cost shift matrix mult
[B96,B97]
0.40
1.62
predictions
115
Cost total tomo
14.00
43.87
116
Cost H mult [B33]
3.02
2.85
E84
117
Prefitting Cost grand total
17.02
46.72
118
Quadratic Incr. WFE from FD3o6 (nm)
with perfect wind
profile
knowledge
18.00
Negative sign indictated reduced
WFE
119
Quadratic Incr. WFE from FD3o6 (nm)
with no wind profile
knowledge
()
60.00
42.00
120
Quadratic Incr. WFE from FD3o6 (nm)
with 20 deg wind
direction
error and 10% wind speed error
41.00
The last column contains the formulas used for DKFo6P (letters
refer to columns and numbers to rows).
Table 8.
DKF Tomography Algorithm Computational Cost Analysis for the Case of
the Kalman Gain
Multiplication
Implemented via SD Convolutions, for the Case
,
,
Denoted, Respectively,
as DKFP-SD and
DKFo6P-SDa
121
DKF Comps with Kalman Gain Mult Implemented via
SD Convolutions
122
DKFP-SD
DKFo6P-SD
DKFo6P-SD Formulas
123
nK
60
120
1D nb of points in convolution
kernels
124
Ng/1e6
0.03
0.12
125
Cost grad est. error
[B90]
3.41
3.16
D106
126
Cost tip/tilt bias removal
[B91]
0.06
C69
127
Cost SD Kalman gain mult
[B93]
669.60
10713.60
128
Cost SD Kalman gain mult/Cost FD
Kalman gain mult
75.68
283.81
D127/D111
129
Mem SD Kalman gain
0.99
3.96
The last column contains the formulas used for DKFo6P-SD (letters
refer to columns and numbers to rows).