Lanczos iteration

Apart from regularization methods, we will also need a method from numerical linear algebra. A comprehensive introduction on this topic can be found, e.g., in [164, 89, 66], on which our discussion in this subsection is based.

The Lanczos iteration is used to find eigenvalues for large matrices, i.e., solving Ax“λx.

This method was introduced in [94, 95] and computes eigenvalues of a real Hermitian matrix quickly if they are well separated from the rest of the spectrum.

For some linear operator A : X Ñ X, and some element x P X, the k-th Krylov subspace K_kpx, Aq is defined as

K_kpx, Aq “tx, Ax, A²x, . . . , A^k´1xu.

The idea of Lanczos iteration is to find a factorization of a matrixAPR^mˆm

A“Q_nT_nQ^T_n, (3.6)

where Q_n P R^mˆn has orthonormal columns and T_n P R^nˆn is a symmetric and Hes-senberg matrix, thus tridiagonal. We denote byα_i the diagonal values of T_n, and the entries of first off-diagonalsβ_i, where iP t1, . . . , nu. In order to obtain a fast method, not the full factorization, i.e.,m “n, is computed, but only a small dimensional sub-space is used, i.e.,n !m.

The factorization (3.6) can be computed by Lanczos iteration, shown in Algorithm 3.1, using an arbitrary starting vectorb. It can be shown, that the columns ofQ_n form an orthonormal basis forK_npb, Aq, with the starting vector b.

42 CHAPTER 3. MATHEMATICAL PRELIMINARIES

Algorithm 3.1 Lanczos Iteration (from [164]) β₀ “0,q₀ “0,b “arbitrary, q₁ “b{}b}.

for n “1,2,3, . . . v “Aq_n α_n“q_n^Tv

v “v´β_n´1q_n´1´α_nq_n β_n “ }v}

q_n`1 “v{β_n

Note that this iterative procedure results in a three-term recurrence only for a Her-mitian (or self-adjoint) matrix, i.e., a matrix that equals its own transposed complex conjugate, for A P C : A^H :“ A¯^T “ A. The three-term recurrence is also due to the orthogonality of the Krylov vectors. Otherwise, in the more general case, the result is ann-step recurrence and the algorithm is then called Arnoldi iteration.

The Lanczos iteration gives a good approximation of eigenvalues of A, which do not cluster, through the eigenvalues of T_n.

When carried out in exact arithmetic, the following theorem, summarizing the prop-erties of the Lanczos iteration, holds:

Theorem 3.1. [164, Theorem 36.1] The matrices Q_n of vectors q_n generated by the Lanczos iteration are reduced QR factors of the Krylov matrix

K_n““

b|Ab| ¨ ¨ ¨ |A^n´1b‰

“Q_nR_n. The tridiagonal matrices T_n are the corresponding projections

T_N “Q^˚_nAQ_n, and the successive iterates are related by the formulae

AQ_n“Q_n`1T˜_n,

where T˜_n is thepn`1q ˆn upper-left section ofT<sub>n`1</sub>, which we can rewrite in the form of a three-term recurrence at step n,

Aq_n“β_n´1q_n´1`α<sub>n</sub>q<sub>n</sub>`β_nq_n`1.

As long as the Lanczos iteration does not break down (i.e., K_n is of full rank n), the characteristic polynomial of T_n is the unique polynomial pⁿ P Pⁿ that solves the Lanczos approximation problem, i.e., that achieves

}pⁿpAqb} “minimum, where Pⁿ “ tmonic polynomials of degree nu.

We will use the Lanczos iteration as basis for a fast blind deconvolution method in Chapter 8. In our application, the main advantage is that invertingA can be reduced to inverting T_n, which is cheap due to the small dimension and the structure of T_n.

CHAPTER 4. A CUMULATIVE RECONSTRUCTOR ON FINITE ELEMENT

BASIS 43

Chapter 4

A cumulative reconstructor on finite element basis

The reconstruction of incoming wavefronts from the measurements made by a WFS is the core of an AO system, performed in the RTC unit. For the RTC of an AO system various methods exist. Most common methods in AO systems operating nowadays are so called matrix-vector-multiplication methods (MVM), which rely on a command matrix mapping WFS measurements directly to actuator commands for the DM. The class of to-be-built telescopes requires new algorithms due to the increase of compu-tational complexity when current methods are scaled to much bigger systems.

Therefore, a new class of algorithms was developed within the last years. These algo-rithms range from direct methods for SCAO to iterative methods for more complex AO systems such as MCAO. Among the direct methods fast algorithms such as CuReD [180, 140] and a hierarchical wavefront reconstruction algorithm [12] were introduced and already tested on sky [13]. For wavefront reconstruction from Shack-Hartmann WFS data several methods are available [170, 124, 65, 37, 115]. As also Pyramid WFS are under consideration for the upcoming ELTs, the number of wavefront reconstruc-tion alogrithms for this sensor type is growing [153, 122, 126].

For more complex AO systems, relying on atmospheric tomography, iterative meth-ods as the Fractal Iterative Method (FrIM) [160], the Finite Element Wavelet Hybrid Algorithm (FEWHA) [179], a gradient-based approach [149, 151] and many others based on Fourier transform or Tikhonov regularization with different penalization terms [177, 53, 64, 132] were developed.

In this chapter, we introduce an alternative version of a cumulative reconstructor acting on a finite element basis (FinECuReD) from [117] and discuss the simulation results obtained with this reconstruction algorithm in ESO’s end-to-end simulation tool OCTOPUS [100], presented in [174], a joint work with Andreas Neubauer and Ronny Ramlau.

44 CHAPTER 4. A CUMULATIVE RECONSTRUCTOR ON FINITE ELEMENT BASIS

4.1 Wavefront sensor geometry

The algorithm is set up for general apertures ΩĂ rA, Bs ˆ rC, Ds without an obstruc-tion. We consider Shack-Hartmann WFS only, thus the aperture Ω is a collection of small subapertures, These subapertures are all of the same quadratic size. Note that for this algorithm no symmetry of Ω is necessary. However, we only consider regions Ω whose boundaries are parallel to the axes. Mathematically this reflects the property that the intersection of each line parallel to an axis is an interval. As an example, we could consider the region shown in Figure 4.1, a quarter of a usual telescope aperture.

Figure 4.1: WFS geometry without obstruction, from [117].

In order to be able to state the algorithm, we need some formalization of the geomet-rical setting above. Let us define a subaperture Ω_ij as follows:

Ω_ij “ rx_i´1, x_is ˆ ry_j´1, y_js, pi, jq PM, where

xi :“A`ih, yj :“C`jh, B :“A`Nxh, D:“C`Nyh,

for some h ą0 and Nx, Ny P N. We denote the set of admissible indices pi, jq byM, i.e.,

M:“tpi, jq PNˆN: 1ďj ďN_y, p

j`1ďiďp_ju

“tpi, jq PNˆN: 1ďiďN_x, q

j`1ďiďq_ju,

CHAPTER 4. A CUMULATIVE RECONSTRUCTOR ON FINITE ELEMENT

BASIS 45

for some 0ďp

j ďp_j ďN_x and 0ďq

j ďq_j ďN_y. Furthermore, let k, l, σ_n, τ_n PN be such that

0“:σ₀ ăσ₁ ă ¨ ¨ ¨ ăσ_k :“N_y, 0 :“τ₀ ăτ₁ ă ¨ ¨ ¨ ăτ_l:“N_x, with

pσn´1 ‰p

σn or p_σn´1 ‰p_σn, 1ănďk, pj “p

σn and p_j “p_σn, σ_n´1 ăj ďσ_n, 1ăn ďk, qτn´1 ‰q

τn or q_τn´1 ‰q_τn, 1ănďl, qi “q

τn and q_i “q_τn, τn´1 ăiďτn, 1ăn ďl.

This description is very technical, but the aperture shown in Figure 4.1 gives a good example, with the following numbers:

Nx “9“Ny, p

1 “5, p

2 “3, p

3 “1, p

4 “p

5 “1, p6 “p

7 “p

8 “p

9 “0, p_i “9,1ďiď7, p₈ “8, p₉ “7, k“7, σ1 “1, σ2 “2, σ3 “3, σ4 “5, σ5 “7, σ6 “8, and, q

i “p

i, q_i “p_i, τ_i “σ_i, l“k, due to symmetry.

Furthermore, we define

L_y :“txPR:px, yq PΩu “ rapyq, bpyqs, yP rC, Ds, L_x :“tyP R:px, yq PΩu “ rcpxq, dpxqs, yP rC, Ds.

Obviously, the functionsa, b, c, d are piecewise constant and it holds that apyq “A`p

σnh, bpyq “ A`p<sub>σn</sub>h, yP pyσn´1, yσnq, 1ďnďk, cpxq “C`q

τnh, dpxq “ C`q_τnh, xP px_τn´1, x_τnq, 1ďnďl.

At jumps, we can only compute one-sided limits, defined by f^`ptq:“ lim

sÑt^`

fpsq and f^´ptq:“ lim

sÑt^´

fpsq.

Im Dokument From Adaptive Optics systems to Point Spread Function Reconstruction and Blind Deconvolution for Extremely Large Telescopes / eingereicht von Dipl.-Ing. Roland Wagner, Bakk.techn. (Seite 53-57)