Optimization of Iterative Control Algorithms for High-Resolution Adaptive Optics Systems

As the number of wavefront sensor subapertures and deformable mirror actuators in adaptive optics systems increases, the computational time of the direct gradient wavefront control algorithm is excessively long, which is a major factor affecting the control performance of adaptive optics systems. The paper combines preprocessing techniques with sparse matrix multiplication techniques to reduce the computational complexity of the algorithm. And the convergence and wavefront control stability of the iterative algorithm are optimized. For an adaptive optical system with 1201 actuators, the computational efficiency of the proposed algorithm is increased by 5 times. And the larger the scale of the adaptive optics system, the more significant the improvement in efficiency compared to the direct gradient wavefront control algorithm.


Optimization of Iterative Control Algorithms for
High-Resolution Adaptive Optics Systems Shenghu Liu , Wang Zhao, Shuai Wang , Kangjian Yang, Ping Yang , Hongli Guan, Han Guo, and Ruifeng He Abstract-As the number of wavefront sensor subapertures and deformable mirror actuators in adaptive optics systems increases, the computational time of the direct gradient wavefront control algorithm is excessively long, which is a major factor affecting the control performance of adaptive optics systems.The paper combines preprocessing techniques with sparse matrix multiplication techniques to reduce the computational complexity of the algorithm.And the convergence and wavefront control stability of the iterative algorithm are optimized.For an adaptive optical system with 1201 actuators, the computational efficiency of the proposed algorithm is increased by 5 times.And the larger the scale of the adaptive optics system, the more significant the improvement in efficiency compared to the direct gradient wavefront control algorithm.

I. INTRODUCTION
A DAPTIVE optics is considered an effective real-time com- pensation technique for mitigating the effects of atmospheric turbulence and thermal blooming, leading to significant improvements in beam quality and enhanced image resolution [1], [2], [3].The direct gradient wavefront control algorithm (DGWC) [4] is a widely used wavefront control algorithm in adaptive optics systems.This algorithm mainly involves the matrix-vector multiplication process.Telescopes have been developing towards large apertures, such as the American 8meter Gemini Telescope [5], the Japanese 8.2-meter Subaru Telescope [6], the American 10-meter Keck Telescope [7], the 30-meter Telescope (TMT) under construction in the United States [8], the 42-meter Extremely Large Telescope (E-ELT) under construction in Europe [9], etc.The adaptive optics system is a necessary configuration of large telescopes.As the scale of adaptive optics systems continues to expand, the number of actuators in the system is showing a rapid growth trend.For example, the PALM-3000 adaptive optics system has more than 3000 actuators of deformable mirrors used to correct high-order aberrations.However, the computational cost of the direct gradient algorithm significantly escalates, with a computational complexity ranging from O(n 2 ) to O(n 3 ) [10].Due to the real-time requirements of adaptive optics systems, the extensive computational load renders the current data processing systems incapable of implementing this algorithm.To address this issue, Luc [11], [12] proposed a multigrid conjugate gradient algorithm for large-scale multi-conjugate systems.However, this method primarily optimizes the computational efficiency of wavefront reconstruction and does not address improvements in the computational efficiency of wavefront control.Chen et al. [10], [13] proposed an iterative control algorithm that uses the sparse characteristics of the slope response matrix, reducing the computational costs of wavefront control.They proved that the conjugate gradient wavefront control algorithm (CGWC) [14], [15] has advantages over other algorithms in terms of computational costs and storage space.But they did not take into account the impact of iterative convergence speed on the computational costs of wavefront control.As the convergence speed of the CGWC algorithm is affected by the spectral properties of the iterative control matrix, slow convergence or non-convergence may occur when the spectral properties of the iterative control matrix are poor [16].Moreover, existing research results do not discuss the stability of iterative wavefront control.This paper focuses on optimizing the iterative convergence speed, computational efficiency, and stability of wavefront control.Based on the CGWC algorithm, the proposed the sparse approximate inverse conjugate gradient wavefront control algorithm (SPAICGWC) combines sparse matrix multiplication techniques [17], [18], [19], [20] with preconditioning techniques [21] to improve computational efficacy.And we ensure the convergence of the iteration and improve the stability of the wavefront control by selecting and modifying the eigenvalues of the iterative control matrix.Simulations and experiments show the effectiveness of the proposed algorithm.For an adaptive optical system with 1201 actuators, the computational efficiency of the proposed algorithm is 6 times that of the DGWC algorithm.The structure of this paper is as follows: Section II presents a detailed introduction to the construction process and principles of the SPAICGWC algorithm, along with an analysis of the algorithm's computational cost.In Section III, sparse threshold and eigenvalue correction value were separately discussed for their impact on the root mean square error of wavefront residual, the number of non-zero elements, and the stability of algorithm iteration counts.Comparative analysis was conducted between the proposed algorithm and the direct gradient wavefront control algorithm as well as the conjugate gradient wavefront control algorithm across nine adaptive optics systems.Section IV validates the algorithm through experimental verification on an adaptive optics system with 169 actuators.Section V concludes with discussions and future prospects.

A. Algorithm Principles and Theoretical Derivation
The DGWC algorithm establishes the relationship matrix between the deformable mirror actuators and the sub-apertures through the influence of each actuator unit voltage on the subaperture slope data.Assuming that a voltage V is applied to the j-th actuator, the average sub-aperture slope of the wavefront can be expressed as: where G x and G y are the average sub-aperture slopes, m is the number of sub-apertures in the wavefront sensor, n is the number of deformable mirror actuators, s i is the normalized area of the sub-aperture, and R j (x, y) is the influence function of the j-th deformable mirror actuator.Because within a certain range, there is a linear relationship between the sub-aperture slope and the actuator voltage, (1) can be represented as: where R xy is the slope response matrix obtained through the deformable mirror and wavefront sensor.In the wavefront control process of an adaptive optics system, control voltages are calculated through the generalized inverse matrix of the slope response matrix R xy , and the control relationship is as follows: where R + xy is the control matrix and is a dense matrix, x is the calculated voltage loaded on the actuators, and g is the detected wavefront slope by the wavefront sensor.By transforming (3), it can be expressed as: The obtained (4) represents the equation for iterative wavefront control, where A is the iterative control matrix and b is the iterative vector.Since the iterative control matrix A is a sparse matrix, the CGWC algorithm calculates the initial search direction p 0 based on (5) with a given initial voltage x 0 : where r 0 is the residual vector.Through the initial search direction and the residual vector, the voltage is iteratively solved, as follows [14]: where P k+1 , r k+1 , and x k+1 are the (k+1)th search direction, residual vector, and updated voltage, respectively.α k and β k are intermediate variables.
The control matrix A obtained from ( 4) is symmetric, but the matrix A is not necessarily positive-definite, so the CGWC algorithm may not necessarily converge.And if matrix A has poor spectral properties, the convergence speed will be slow.
Therefore, the CGWC algorithm needs to be optimized.The main steps of the proposed algorithm are as follows: Step 1: Transformation of the linear system using preconditioning techniques.
Step 2: Sparsification of the control matrix.
Step 3: The selection and correction of eigenvalues of the sparsified control matrix.
Step 4: Iterative solution of the final generated linear equations using the conjugate gradient algorithm.
First, the linear system is transformed.For (4), the following equation can be constructed: If a preprocessing matrix D can be found that A ai possesses better spectral properties and can be sparsified, then combining the advantages of preprocessing and sparse matrix multiplication can enhance the convergence speed of the iterative wavefront control algorithm.Currently, methods commonly used for constructing preprocessing matrices include banded preconditioning, triangular preconditioning, incomplete decomposition preconditioning, and sparse approximate inverse [21], [22], [23], [24].Sparse approximate inverse [25], [26], [27] is a method to construct an approximate inverse matrix with a sparse structure based on the sparse characteristics of a matrix.Therefore, the preprocessing matrix D can be constructed using sparse approximate inverses.For the sparse approximate inverse algorithm, matrices D and A satisfy the following relationship: problems.The d j and e j are the column vectors of matrices D and I, respectively.Then, a sparse structure is defined for matrix D, denoted as matrix Z, where Z elements consist of 0 and 1, with 1 representing non-zero elements and 0 representing zero elements.In actual calculations, (8) can be rewritten as where is the residual threshold.When ( 9) is satisfied, A d j = e j is considered.Assuming an initial diagonal matrix Z as the sparse structure matrix, in order to satisfy the residual boundary conditions of ( 9), non-zero positions can be added to d j based on the sparse structure matrix of matrix A. Furthermore, from ( 7) and ( 9), it can be deduced that the preprocessing matrix D is the approximate inverse matrix of matrix A. The matrix A ai is an approximate unit matrix of I, with non-zero elements mainly concentrated on the main diagonal of A ai and a relatively large proportion of approximate zero elements.Therefore, matrix A ai can be further sparsified, utilizing sparse matrix multiplication to enhance computational efficiency.By sparsifying matrix A ai through a sparse structure matrix W, the resulting new sparse matrix A spai can be represented as: (10) where represents dot multiplication, and W is a matrix with a 0 and 1 distribution.
Next, the matrix A spai is an approximate unit matrix of I, with some eigenvalues fluctuating around 1. Therefore, the spectral properties of A spai are better than those of matrix A. However, there may be several eigenvalues that are approximately 0 or significantly less than 1.Correcting these eigenvalues can ensure convergence and improve the convergence speed.The SVD decomposition is performed on the matrix A spai as follows: After sorting the diagonal matrix S in descending order, eigenvalue selection and correction are performed as follows: where ω is the correction value, S (i,i) is the diagonal element of matrix S, and S T is the screening threshold.S T is determined using the largest gradient of adjacent eigenvalues, as follows: Substituting the selected and corrected eigenvalues into (11), the sparse approximate inverse wavefront control matrix is represented by (14).This ensures the sparsity of the preprocessing matrix D and the favorable spectral properties of the new matrix A new , effectively enhancing the convergence speed of the algorithm.
By substituting ( 14) into (7), the iterative solution of the sparse approximate inverse wavefront control algorithm can be expressed as follows: where y is the intermediate variable.
Finally, the voltage is solved utilizing the conjugate gradient algorithm as in iteration process (6).
Combining the derived formulas, the overall process of designing the algorithm is as follows: Firstly, construct the preconditioning matrix D using (8), and use the preconditioning matrix to transform the linear system, obtaining the iterative control matrix Secondly, apply sparsification to the iterative control matrix A ai using (10) to obtain the matrix Then, modify the eigenvalues of the matrix A spai using (12).Finally, solve for the voltage x through the iterative process described in (6).

B. Computational Cost of Algorithms
The computational process of the DGWC algorithm is described by (3), where the control matrix R + xy is a dense matrix with n rows and 2m columns.Therefore, the computational cost for solving the voltage x through matrix-vector multiplication is as follows: The CGWC algorithm iteratively solves the voltage x , based on (4) and the calculation process is shown in (6).When x k+1 ≈ x k , then x k+1 is the solution of (4).Since the iteration control matrix A and slope response matrix R are sparse matrices, the dominant part of the calculation cost of sparse matrix multiplication is calculated as follows: where k represents the iteration count, and c is the number of non-zero elements in the iteration control matrix A.
Based on (4), the proposed algorithm, during the execution of linear system transformation, sparsification of the iterative control matrix, and eigenvalue correction, introduces modifications to the iterative control matrix.Additionally, it introduces extra sparse matrix-vector multiplication, as indicated by (15).The proposed algorithm solves the voltage through (15), and the solving process is consistent with (6) in the CGWC algorithm.The computational cost of this multiplication is as follows: where c 1 represents the number of non-zero elements in the new iteration control matrix A new , and c 2 ; represents the number of non-zero elements in the preprocessing matrix D.

III. COMPARISON OF ALGORITHM PERFORMENCE
Table I presents data for nine adaptive optics systems, with specific parameters.The c represents the number of non-zero elements in the iterative control matrix of the CGWC algorithm.The computational cost of the algorithm can be calculated using Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.16), (17), and (18).And the stability of the algorithm's iteration count is provided by (19):

TABLE I PARAMETERS OF NINE ADAPTIVE OPTICS SYSTEMS
where k represents the average iteration number, T signifies the number of closed-loop control cycles in the adaptive optics system, and ε denotes the mean square deviation of the algorithm's iteration count.
Taking the adaptive optics system with a deformable mirror of 265 actuators as an example, the distribution of nonzero elements in the preprocessing matrix D is shown in Fig. 1(a).It can be observed that matrix D is sparse.The eigenvalue distribution of the control matrix A ai obtained through the preprocessing matrix D is shown in Fig. 1(b), indicating the presence of a few eigenvalues significantly smaller than 1.The sparse structure and distribution of antidiagonal elements of matrix A ai are illustrated in Fig. 2. Despite a considerable proportion of nonzero elements in the matrix A ai , as illustrated in Fig. 2(a), the distribution of its anti-diagonal elements, as shown in Fig. 2(b), suggests a substantial number of elements approximating zero.This observation aligns with the analysis in Section II.

A. The Impact of Sparsity Threshold
From ( 18), it can be deduced that the number of iterations in the algorithm and the number of nonzero elements in the control matrix A new are the primary factors influencing the computational cost of the SPAICGWC algorithm.Therefore, an analysis was conducted to examine the impact of the sparsity level of the control matrix A new on the SPAICGWC algorithm on the adaptive optics system with 256 actuators.
Fig. 3(a) illustrates the trends in wavefront residual RMS and the number of nonzero elements as the sparsity threshold increases.When the sparsity threshold ranges from 0 to 0.008, the wavefront residual RMS exhibits a decreasing trend, while in the range of 0.008 to 0.02, the wavefront residual shows an increasing trend.The number of nonzero elements gradually decreases as the sparsity threshold increases, with a more pronounced decrease when the sparsity threshold is between 0 and 0.01.Beyond a sparsity threshold of 0.01, the decrease in the number of nonzero elements becomes more gradual.Fig. 3(b) presents the trends in wavefront residual RMS and algorithm iteration stability as the sparsity threshold increases.The algorithm iteration stability gradually improves with an increase in the sparsity threshold, particularly noticeable in the range of sparsity thresholds from 0 to 0.01.When the sparsity threshold exceeds 0.01, the improvement in algorithm iteration stability becomes less pronounced.This indicates that selecting an appropriate sparsity threshold not only ensures the accuracy of wavefront control but also effectively reduces the computational cost of the algorithm while enhancing the stability of algorithm iterations.

B. The Impact of Eigenvalue Correction Values
Eigenvalue correction involves further adjustments to the control matrix after sparsification, primarily affecting the wavefront residual RMS and the stability of algorithm iteration.The control matrix after sparsification for the adaptive optics system with 256 actuators was subjected to eigenvalue correction.And 400 different correction values are tested, with each correction value subjected to 100 control experiments in each group, as depicted in Fig. 4. Wavefront residual RMS initially decreases and then Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.increases with an increase in the eigenvalue correction value.Specifically, when the eigenvalue correction value is in the range of 0 to 0.2, the wavefront residual exhibits a decreasing trend, while in the range of 0.2 to 0.4, it shows an increasing trend.
Furthermore, the algorithm iteration stability gradually improves when the eigenvalue correction value falls within the ranges of 0 to 0.22 and 0.3 to 0.4.However, the algorithm iteration stability deteriorates when the eigenvalue correction value is in the range of 0.22 to 0.3.
As a result, there exist specific eigenvalue correction values that lead to minimal values for both wavefront residual RMS and algorithm iteration stability.This indicates that by carefully balancing the wavefront residual and algorithm stability and selecting appropriate eigenvalue correction values for the modified control matrix after sparsification, further improvements can be achieved in wavefront control accuracy and algorithm iteration stability.

C. Simulation Validation
Testing was conducted on the remaining 8 sets of adaptive optical systems listed in Table II.Generate 800 sets of aberrated wavefronts, with 400 used to choose sparse thresholds and eigenvalue correction values for the proposed method and another 400 used to compare algorithm performance.The generated distorted wavefront conforms to the Kolmogorov turbulence model and satisfies d/ R 0 = 1 [28], where d is the wavefront sensor subaperture diameter and R 0 is the atmospheric coherence length.The selection of the sparsity threshold and eigenvalue correction values is consistent with the approach used for the adaptive optics system with 256 actuators.Table II presents the parameter selection and the corresponding number of nonzero elements in the control matrix and the approximate inverse matrix using the SPAICGWC algorithm for nine sets of adaptive optics systems.To provide a visual representation of the algorithm's computational cost, as shown in Fig. 5, when the number of actuators of adaptive optics system exceeds 265 elements, the computational cost of the CGWC algorithm is lower than that of the DGWC algorithm.In the case of the proposed algorithm, the computational cost is lower than the DGWC algorithm when the number of deformable mirror actuators exceeds 61 elements.Fig. 5 indicates that for adaptive optical systems with actuators ranging from 61 to 1201 elements, the computational cost of the SPAICGWC algorithm is consistently lower than that of the CGWC algorithm.Additionally, for the adaptive optical system with 1201 actuators, the proposed SPAICGWC algorithm offers a sixfold improvement in computational efficiency compared to the DGWC algorithm.Therefore, the computational efficiency of the proposed algorithm is higher than that of the CGWC and DGWC algorithms.There are two main reasons for the improvement of the computational efficiency of this algorithm.One using the preprocessing methods can speed up the iterative convergence speed and reduce the number of iterations to solve the voltage, as shown in Fig. 6(a).The other, the constructed preprocessing matrix is a sparse matrix, so sparse matrix multiplication can be used during the calculation to reduce the computational cost.The advantage of the SPAICGWC algorithm in terms of computational cost becomes more significant with an increasing number of actuators in the adaptive optics system.Fig. 6(a) shows the convergence speed of the CGWC and the SPAICGWC algorithms, and the smaller the number of iterations, the faster the convergence speed of the algorithm.And the iteration count of the SPAICGWC algorithm is significantly lower than that of the CGWC algorithm.As the actuator scale and Shack-Hartmann spatial resolution of the adaptive optics system increase, the number of iterations of the SPAICGWC algorithm gradually increases.However, the number of iterations using the CGWC algorithm does not gradually increase with the number of actuators, as observed in adaptive optics systems with deformable mirrors with 61, 109, 321, and 713 actuators.This is because the unfavorable spectral properties of the slope response matrix slow down the convergence speed of the CGWC algorithm.The SPAICGWC algorithm provides the linear system with favorable spectral properties and the same solution.Correcting the eigenvalues can further improve the spectral properties of the matrix.The proposed algorithm effectively avoids this issue, so the convergence speed of the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.proposed algorithm mainly depends on the number of actuators in the adaptive optics system.Furthermore, when the number of adaptive optical system actuators is less than 321, the proposed SPAICGWC algorithm exhibits iteration stability consistently less than 1, whereas the iteration stability for the CGWC algorithm is consistently greater than 1, as depicted in Fig. 6(b).The improvement in the wavefront control stability of the proposed algorithm is mainly due to the sparseness of the iterative control matrix.Matrix sparsification can effectively suppress the impact of system noise on iterative control.
For the wavefront control experiments conducted on the nine sets of adaptive optics systems, the proposed algorithm exhibits control accuracy comparable to the DGWC and CGWC algorithms.The average wavefront residual RMS for each set of adaptive optics systems is shown in Fig. 7(a).The control accuracy of the proposed algorithm is the same as that of the CGWC and DGWC algorithms.Fig. 7(b) shows that the SPAICGWC control process is consistent with CGWC and DGWC.Simulation results demonstrate that the proposed SPAICGWC algorithm improves computational efficiency and algorithm iteration stability while ensuring control accuracy consistent with DGWC and CGWC under the same conditions.

IV. EXPEREMENT
We conducted experimental verification on an adaptive optics system with 169 elements.Fig. 8(a) depicts the experimental optical setup, where collimated light passes through a tilt mirror and is reflected by a deformable mirror.The specific method of the experiment is as follows.First, the slope response matrix R xy is obtained by applying a unit voltage to each deformable mirror actuator.Then, two sets of different slope data were obtained using two aberration plates; one set was used to select the sparse threshold and eigenvalue correction values, and the other set was used for verification.Fig. 9(a) illustrates the distribution of the sparse structure of the slope response matrix for the 169-actuator adaptive optics system.The eigenvalues and their selection distribution are presented in Fig. 9(b).In the experiments, the voltages obtained by solving the DGWC algorithm were used as the estimated true values to determine the sparse threshold and eigenvalue correction values for the SPAICGWC algorithm.In Fig. 10, the curve depicting voltage solution accuracy and the variation of nonzero elements with sparse threshold indicates an optimal sparse threshold of 0.009, with an optimal eigenvalue correction value of 0.375.
The experimental results are shown in Table III.The computational efficiency of the SPAICGWC algorithm is twice that of the DGWC algorithm and three times that of the CGWC algorithm.In the experimental results, the non-zero elements of the iterative control matrix of the CGWC algorithm range between 109 and 265 executors in Table I.The non-zero elements of the iterative control matrices A new and matrix D for the SPAICGWC algorithm are between 109 and 265 executors in Table II.This shows that there is consistency between the experimental and simulated results.
Fig. 11(a) depicts the voltages applied to the deformable mirror actuators obtained from the wavefront control algorithms.The results from the CGWC algorithm, the SPAICGWC algorithm, and the DGWC algorithm exhibit consistent outcomes.Using the three algorithms for wavefront control, the control execution time is presented in Fig. 11(b).Notably, the execution time of the SPAICGWC algorithm is shorter compared to the DGWC and the CGWC algorithms.

V. CONCLUSION
In this paper, we have introduced a sparse approximate inverse wavefront control algorithm.Numerical simulations have demonstrated that as the scale of the adaptive optics system increases, the computational cost of the proposed algorithm is significantly lower than that of the direct gradient wavefront control algorithm and the conjugate gradient wavefront control algorithm.Regarding solution stability, the iterative stability of the sparse approximate inverse conjugate gradient wavefront control algorithm superior to that of the conjugate gradient wavefront control algorithm.
Furthermore, through validation on an adaptive optics system with a deformable mirror consisting of 169 actuators, the effectiveness of the algorithm has been confirmed, showing a twofold improvement in computational efficiency compared to the direct gradient wavefront control algorithm.Therefore, the proposed algorithm is expected to offer a feasible solution for accelerating computations in adaptive optics systems with more than a thousand actuators.

Fig. 1 .
Fig. 1.Preprocessing matrix and generated control matrix for 265 actuators.(a) The nonzero element distribution of the preprocessing matrix D. (b) The eigenvalue distribution of the control matrix.

Fig. 2 .
Fig. 2. Non-zero element and anti-diagonal element distribution of the control matrix A ai for 265 actuators.(a) Nonzero element distribution.(b) The antidiagonal element distribution.

Fig. 3 .
Fig. 3.The distribution of wavefront residual RMS, nonzero elements, and stability of algorithm iteration with respect to the sparse threshold value δ.(a) Distribution of wavefront residual RMS and nonzero elements.(b) Distribution of wavefront residual RMS and stability of algorithm iteration.

Fig. 4 .
Fig. 4. Impact of eigenvalue correction values on wavefront residual RMS and stability of iteration counts.

Fig. 6 .
Fig. 6.Comparison of iteration counts and iteration stability for iterative control algorithms.(a) Comparison of average iteration counts between CGWC and SPAICGWC algorithms.(b) Comparison of iteration count stability between CGWC and SPAICGWC algorithms.

Fig. 7 .
Fig. 7. Comparison of wavefront residual RMS among nine adaptive optics systems.(a) Comparison of average wavefront residual RMS.(b) Comparison of wavefront residual RMS during the control process.

Fig. 9 .
Fig. 9. Non-zero element distribution of control matrix A ai along with the generated control matrix's eigenvalue distribution.(a) Distribution of non-zero element.(b) Distribution and selection of eigenvalues.

Fig. 10 .
Fig. 10.The impact of sparsity threshold and eigenvalue correction on voltage solution accuracy, control matrix nonzero elements, and algorithm iteration stability.(a) The influence of sparsity threshold on voltage solution accuracy and nonzero element count.(b) The impact of eigenvector correction values on voltage solution accuracy and iteration stability.

TABLE II PARAMETERS
OF NINE SETS OF ADAPTIVE OPTICS SYSTEMS FOR SPARSE APPROXIMATE INVERSE WAVEFRONT CONTROL ALGORITHM

TABLE III ALGORITHM
PERFORMANCE COMPARISONAuthorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.