1.

Introduction

Future space coronagraphs attempting to image and characterize Earth-like planets around nearby solar-type stars will rely critically on closed-loop wavefront sensing and control (WFS&C) using deformable mirrors (DMs) to mitigate optical aberrations in the primary beam path. These aberrations, primarily mid-spatial frequency wavefront errors and optical misalignments in the telescope and coronagraph optics, give rise to a speckle floor that is coherent with the star and evolves slowly over time in response to miniscule drifts in the thermal and mechanical state of the observatory. If uncorrected, the speckle floor overwhelms the faint image of the orbiting planet, which is expected to be ${10}^{10}$ times fainter than the host star at 0.1 arcseconds of separation or less at visible wavelengths.¹

The current state-of-the-art wavefront control algorithms, stroke minimization (SM)² and electric field conjugation (EFC),³ compute DM command solutions using a first-order approximation of the focal-plane electric field in an optimal control framework. In both cases, the optimal DM update is written down in closed form as the solution of a linear system of equations constructed from a Jacobian matrix that describes the influence of each DM actuator on the focal-plane electric field.

Until recently, little attention has been paid to the computational demands of SM and EFC. The complexity of SM and EFC is dominated by the cost of computing and manipulating the Jacobian matrix, whose size is proportional to the product of the number of DM actuators and number of dark-zone pixels. In broadband imaging, this is compounded by the requirement to compute a separate Jacobian for each controlled wavelength. The Jacobian is most often model-based, in which case an optical diffraction model of the coronagraph is evaluated repeatedly to predict the focal-plane influence of each actuator. Moreover, the Jacobian is a linearization of the true, nonlinear behavior of the DMs and must be recalculated periodically as the state of the DMs evolves over time.

As direct imaging missions demand DMs with ever-higher actuator density to enable wider and wider search areas, computational aspects will inevitably become a point of concern from a systems engineering standpoint. In on-orbit WFS&C, all sensing and control computations are processed by the flight computer; conversely, in ground-in-the-loop scenarios, raw data is communicated to a ground-based computing node that calculates the DM correction and relays the correction back to the observatory. Though each approach has tradeoffs, a major advantage of on-orbit WFS&C is the ability to update DM commands more frequently without relying on the continuous availability of a communication link with the ground station. On-orbit WFS&C can help to relax observatory-level wavefront stability requirements by enabling the high-contrast dark zone to be maintained over shorter time intervals. Successful deployment of an on-orbit architecture is predicated on the availability of sufficient computational resources; however, radiation-hard, space-qualified computing hardware lags behind conventional hardware by decades and is extremely resource-limited. This poses a substantial computation capability gap if the current algorithms are expected to be deployed on-orbit on a future direct imaging mission, such as NASA’s planned Habitable Worlds Observatory (HWO).

As a case in point, the current state-of-the-art testbed for coronagraph laboratory demonstrations, the decadal survey testbed (DST) at NASA’s Jet Propulsion Laboratory, has successfully demonstrated $3.82\times {10}^{-10}$ instrumental contrast over a 10% fractional bandpass within an annular dark zone extending from $3{\lambda }_0/D$ to $8{\lambda }_0/D$, where ${\lambda }_0$ is the central imaging wavelength and $D$ is the diameter of the telescope aperture, using two DMs each with $48\times 48$ actuators.⁴ The baseline requirements for an HWO-like mission will be considerably steeper, using the Large UV/Optical/Infrared Surveyor (LUVOIR) and Habitable Exoplanet Observatory (HabEx) flagship mission concepts formulated for the Astro2020 Decadal Survey^5,6 as representative reference designs. The HabEx design features two $64\times 64$ DMs, and a search area with a maximal outer radius of $32{\lambda }_0/D$ over a 20% bandpass—nearly double the total actuator count, a factor of 4 increase in search radius, and a factor of two increase in control bandwidth. Meanwhile, the LUVOIR Architecture “A” reference design includes a pair of $128\times 128$ DMs and a dark zone with a $64{\lambda }_0/D$ maximal outer radius over a 10% fractional bandpass or wider. Because the number of detector pixels in the dark zone scales with the area of the dark zone rather than its radius, these parameters correspond to an increase by a factor of 32 (Habex) and 512 (LUVOIR-A) of the worst-case Jacobian dimensionality, for each control wavelength, compared to DST.

In recent work, we formulated an alternative wavefront control framework that iteratively compute DM updates using gradient-based optimization techniques, which eliminates the need to calculate and manipulate the Jacobian.⁷ Our approach is based in part on a numerical technique known as algorithmic differentiation (AD),⁸ to calculate the gradients of the cost function for optimal control accurately and efficiently. We described an AD-based counterpart to SM, which we termed AD penalty stroke minimization (AD-PSM), and compared it to SM in simulations with a small-angle LUVOIR design. Our results indicated superior computational efficiency with AD-PSM and comparable starlight suppression performance. While the CPU time and memory consumption of SM grew superlinearly with actuator count, the increase in both for AD-PSM was negligible (e.g., with $128\times 128$ actuators, AD-PSM utilized 95% less memory and CPU time), suggesting that iterative methods are a promising alternative to Jacobian-based techniques for on-orbit WFS&C with high actuator counts and large dark zones.

In this paper, we report on the first experimental demonstration of AD-PSM as well as an AD-based counterpart to the EFC algorithm, termed AD-EFC, using the High-contrast Imager for Complex Aperture Telescopes (HiCAT) at the Space Telescope Science Institute in Baltimore, Maryland. We benchmark the contrast performance of AD-PSM and AD-EFC as a function of several key parameters including regularization and the termination tolerance of the nonlinear optimizer, and compare it to SM and EFC.

This paper is structured as follows. In Sec. 2, we review concepts from our earlier work, including AD and the mathematical formulation of AD-PSM and AD-EFC. In Sec. 3, we provide an overview of HiCAT and discuss our experimental setup, including algorithm implementation details that are pertinent to our demonstration. In Sec. 4, we present and discuss our experimental results. Finally, in Sec. 5, we provide our conclusions and briefly describe our planned future work.

2.

Wavefront Control Using Algorithmic Differentiation

The goal of the WFS&C loop in coronagraphy is to drive starlight within the dark zone toward zero over time so that a faint orbiting companion becomes detectable against the reduced starlight background. Each iteration of the WFS&C loop, indexed by the integer $k$, consists of two steps: (1) an estimation step, in which an estimate ${\widehat{\mathbf{E}}}_{\mathrm{ab},k}$ of the true aberrated electric field ${\mathbf{E}}_{\mathrm{ab},k}$ is formed within the dark zone from focal-plane intensity measurements, and (2) a control step, in which the DM correction is updated to reduce the energy in ${\mathbf{E}}_{\mathrm{ab},k}$, as shown in Fig. 1. In this paper, we focus principally on the control step.

Fig. 1

Closed-loop coronagraphic WFS&C uses images from the science camera, rather than an external wavefront sensing instrument, to estimate the electric field from the host star within the dark zone and drive it toward zero. Model-based WFS&C algorithms use a numerical model of the coronagraph to solve an inverse problem for the unknown electric field and corresponding DM correction, respectively.

Modern model-based wavefront control algorithms find the DM correction update ${\mathbf{a}}_k$ by minimizing some cost function $J_k({\mathbf{a}}_k;{\widehat{\mathbf{E}}}_{\mathrm{ab},k})$ with respect to ${\mathbf{a}}_k$. Usually, $J_k$ is constructed to trade off between minimizing starlight and minimizing the size of the correction, which helps to regularize the problem and stabilize the solution. In general, the true relationship between the starlight in the dark zone and the DM correction is highly nonlinear and nonconvex, owing to the fact that the DMs impart phase-only corrections of the form $\exp \{i{\mathit{\varphi }}_{\mathrm{DM}}\}$ at or near the coronagraph entrance pupil. However, when the optical aberrations are small, we can approximate the true electric field in the coronagraph entrance pupil with a first-order Taylor series expansion about a small update to the DM commands. In this case, the corrected electric field in the dark zone has the form

This is a quadratic function of ${\mathbf{a}}_k$, meaning that under this approximation, there exists a unique, optimal DM correction that minimizes the dark-zone starlight.

SM and EFC utilize the relationship in Eq. (5) to derive closed-form expressions for this optimal correction in terms of ${\mathbf{G}}_k$ that can be evaluated by solving a linear system of dimension $N_{\mathrm{act}}\times N_{\mathrm{act}}$, as illustrated in Fig. 2. As we discuss in Appendix B, this is equivalent to minimizing the cost function using Newton’s method with an exact Hessian matrix. Alternatively, it is possible to find the DM correction by minimizing the cost function $J_k$ with respect to ${\mathbf{a}}_k$ iteratively, rather than analytically, using gradient-based nonlinear optimization as shown in Fig. 3. To do so eliminates the need to calculate the Jacobian matrix, but requires a way of calculating the gradient vector $\partial J_k/\partial {\mathbf{a}}_k$. Reverse-mode algorithmic differentiation (RMAD) provides a way of doing so that is both computationally efficient and accurate, in the sense that the derivatives computed by RMAD are accurate to machine precision and do not utilize finite-difference approximations.⁸ The basic principle of RMAD is that any function that can be written down as a sequence of differentiable operations, called the forward model, can be transformed mechanistically to construct a related function, called the adjoint model, that evaluates the derivative of the forward model with respect to any of the intermediate variables encountered during its evaluation, as well as its inputs. Figure 4 illustrates this procedure for the wavefront control cost function $J_k$. We refer the reader to our earlier work for a more comprehensive discussion of the principles of RMAD and its application to wavefront control.⁷

Fig. 2

In EFC and SM, the Jacobian matrix is precomputed using a computer model of the coronagraph to predict the effect of an update to each of the $N_{\mathrm{act}}$ DM actuators individually on the dark-zone electric field, here represented by the Kronecker $\delta$ functions ${\delta }_n$, where ${\delta }_n[i]=1$ if $i=n$ and zero otherwise. The Jacobian ${\mathbf{G}}_k$ and the estimate of the aberrated dark-zone electric field ${\widehat{\mathbf{E}}}_{\mathrm{ab},k}$ together are used to construct a linear system of equations whose solution is the desired DM update ${\mathbf{a}}_k^{*}$. We show in Appendix B that this is equivalent to minimizing the wavefront control cost function using Newton’s method.

Fig. 3

Our AD-based wavefront controllers AD-PSM and AD-EFC use RMAD to differentiate the wavefront control cost function $J_k$ with respect to the DM correction update ${\mathbf{a}}_k$, yielding the gradient vector $\partial J_k/\partial {\mathbf{a}}_k$ evaluated at the current iterate ${\mathbf{a}}_k^n$. A nonlinear optimization algorithm calculates a new iterate ${\mathbf{a}}_k^{n+1}$ that reduces the value of the cost function, i.e., $J_k({\mathbf{a}}_k^{n+1})\le J_k({\mathbf{a}}_k^n)$. This procedure is repeated until the gradient becomes sufficiently small, indicating that the solution ${\mathbf{a}}_k^{*}$ is at or near a local minimum of the cost function. A starting guess for the solution ${\mathbf{a}}_k^0$ as well as the aberrated electric field ${\widehat{\mathbf{E}}}_{\mathrm{ab},k}$ are the input parameters.

Fig. 4

The forward model for the wavefront control problem maps DM command updates ${\mathbf{a}}_k$ to a scalar cost function $J_k$. The DM command vector is split into independent command vectors ${\mathbf{a}}_{1,k}$ and ${\mathbf{a}}_{2,k}$ for the pupil-plane and out-of-pupil DM, respectively. These are mapped onto DM surfaces ${\mathbf{s}}_{1,k}$ and ${\mathbf{s}}_{2,k}$ using the model in Appendix C, and propagated through an end-to-end coronagraph model to predict the resulting change in dark-zone electric field ${\mathbf{E}}_{\mathrm{DM},k}$. Reverse-mode AD transforms the forward model into an adjoint model that backpropagates the partial derivatives of $J_k$ with respect to each intermediate variable ${\mathbf{a}}_{1,k}$, ${\mathbf{a}}_{2,k}$, ${\mathbf{s}}_{1,k}$, ${\mathbf{s}}_{2,k}$, and ${\mathbf{E}}_{\mathrm{DM},k}$ in reverse order, starting from the output on the right. The derivatives with respect to the individual DM command vectors ${\mathbf{a}}_{1,k}$ and ${\mathbf{a}}_{2,k}$ are concatenated to form the full gradient vector for optimization.

3.

Experimental Setup

In this section, we provide an overview of the HiCAT testbed and provide details about the implementation of AD-PSM and AD-EFC, the reference Jacobian-based implementations of SM and EFC, and the electric field estimation algorithm used in the estimation step of the WFS&C loop.

3.1.

HiCAT Testbed

HiCAT is a testbed dedicated to technology demonstrations for coronagraphy on segmented-aperture space observatories, with the intent of being directly traceable to a future HWO-like mission.^10–17 These technologies include Lyot coronagraphy, high-order WFS&C for generating and stabilizing dark zones, and low-order wavefront sensing (LOWFS).¹⁸ HiCAT operates in a mid-contrast regime (${10}^{-7}$ to ${10}^{-8}$) which approaches the limit achievable outside of a vacuum environment. The testbed is equipped with two Boston Micromachines Kilo-DMs for high-order sensing and control, with 952 actuators each, making it suitable for our proof-of-concept demonstrations. One DM is placed in a plane conjugate to the HiCAT entrance pupil, while the second DM is placed $\sim 300\mathrm{mm}$ farther along the optical axis, corresponding to a Fresnel number $N_F\approx 98$ at a wavelength of 638 nm. This configuration enables simultaneous control of amplitude and phase aberrations over a dark zone that extends over both halves of the image plane. We conducted our experiments using a Thorlabs MCLS1 laser diode source, which emits monochromatic light with a central wavelength ${\lambda }_0=638\mathrm{nm}$.

HiCAT additionally has an IrisAO segmented DM with 37 hexagonal segments, each with controllable piston/tip/tilt, to act as a telescope pupil simulator and to inform experimental efforts devoted to segment-level tolerancing and stabilization.^19,20 Figure 5 shows a simplified system layout of HiCAT, including the primary imaging beam path as well as several additional beam paths used by the LOWFS and metrology subsystems.

Fig. 5

Simplified, partially-unfolded layout of the HiCAT testbed. The elements encountered by the primary imaging beam path are highlighted in bold. DM1 and DM2 indicate the in-pupil and out-of-pupil DMs, respectively. Our experiments utilized a coronagraph configuration in which the apodizer is replaced with a flat mirror (see Sec. 3.1).

Our experiments on HiCAT utilized a classical Lyot coronagraph (CLC) design with a hexagonal entrance pupil mask designed to mask the extreme edges of the IrisAO, a circular focal-plane mask, and a circular Lyot stop. Figure 6 shows an overlay of the CLC pupil masks along with a simulated stellar point-spread functions (PSFs). Figure 7 shows example experimental PSFs obtained before and after closed-loop WFS&C using SM, along with the corresponding DM commands.

Fig. 6

(a) Overlay of HiCAT pupil masks projected onto the in-pupil DM plane, including the reflective area of the DM, the IrisAO segmented aperture, the entrance pupil mask (with the same geometry, but slightly undersized relative to the IrisAO) and the Lyot stop. (b) Simulated and (c) experimental coronagraphic PSFs with outline of the geometrical edge of the focal-plane mask.

Fig. 7

Experimental on-axis images from HiCAT. (a) Non-coronagraphic image, (b) coronagraphic image prior to the WFS&C loop, (c) coronagraphic image after 80 iterations of SM, and corresponding actuator commands for the in-pupil DM (d), and out-of-pupil DM (e). In panels (b) and (c), the geometrical edge of the focal-plane mask (FPM) with radius $3.34{\lambda }_0/D_{\mathrm{LS}}$ is shown as well as the inner and outer edges of the dark zone at $5.8{\lambda }_0/D_{\mathrm{LS}}$ to $9.8{\lambda }_0/D_{\mathrm{LS}}$, respectively, where $D_{\mathrm{LS}}$ is the Lyot stop diameter and ${\lambda }_0=638\mathrm{nm}$.

3.3.

Estimation Algorithm

We used the pairwise probe estimator³ for the estimation step in all experiments. The pairwise estimator forms a least-squares estimate ${\widehat{\mathbf{E}}}_{\mathrm{ab},k}$ of the focal-plane electric field ${\mathbf{E}}_{\mathrm{ab},k}$ by applying a series of $P$ probing DM commands ${\mathbf{u}}_p$ to generate probing electric fields ${\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_p)$ that interfere with ${\mathbf{E}}_{\mathrm{ab},k}$. The data vector for the least-squares estimate is formed by differencing the images resulting from ${\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_p)$ and ${\mathbf{E}}_{\mathrm{DM},k}(-{\mathbf{u}}_p)$. The estimate of the $m$’th pixel in the dark zone is then found by finding the least-squares solution of the system

Eq. (13)

$$\left[\begin{array}{c}{\mathbf{I}}_{k,1}^{+}[m]-{\mathbf{I}}_{k,1}^{-}[m]\\ ⋮\\ {\mathbf{I}}_{k,P}^{+}[m]-{\mathbf{I}}_{k,P}^{-}[m]\end{array}\right]=\left[\begin{array}{cc}\text{Re}\{{\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_1)\}[m]& \mathrm{Im}\{{\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_1)\}[m]\\ ⋮& ⋮\\ \mathrm{Re}\{{\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_P)\}[m]& \mathrm{Im}\{{\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_P)\}[m]\end{array}\right]\left[\begin{array}{c}\mathrm{Re}\{{\mathbf{E}}_{\mathrm{ab},k}\}[m]\\ \mathrm{Im}\{{\mathbf{E}}_{\mathrm{ab},k}\}[m]\end{array}\right],$$

where ${\mathbf{I}}_{k,p}^{\pm }\triangleq {|{\mathbf{E}}_{\mathrm{ab},k}+{\mathbf{E}}_{\mathrm{DM},k}(\pm {\mathbf{u}}_p)|}^2$.

We generated four DM probe functions ${\mathbf{u}}_p$ that were optimized to produce probing fields of the form

Eq. (14)

$${\mathbf{E}}_{\mathrm{DM},k}^p=\mathrm{sgn}\{{\mathit{\theta }}_x\}\exp \{i\pi \frac{p-1}{4}\mathrm{sgn}\{{\mathit{\theta }}_x\}\},$$

where $p\in \{\mathrm{1,2},\mathrm{3,4}\}$ and $\mathrm{sgn}\{{\mathit{\theta }}_x\}$ is the sign of the focal-plane $x$ coordinate. These probing fields, which have anti-Hermitian symmetry, correspond to purely real commands applied to the in-pupil DM, which we can calculate via a regularized least-squares approach for each $p$,

Eq. (15)

$$\underset{{\mathbf{u}}_k^P}{\arg \min }{\Vert {\mathbf{G}}_{1,k}{\mathbf{u}}_k^p-{\mathbf{E}}_{\mathrm{DM},k}^p\Vert }^2+{\alpha }_{\text{probe}}^2{\Vert {\mathbf{u}}_k^p\Vert }^2,$$

where ${\alpha }_{\text{probe}}^2$ is the Tikhonov regularization term for the probe calculation and ${\mathbf{G}}_{1,k}$ is the Jacobian for the in-pupil DM. This is identical to the EFC problem in Eq. (10) but with the probing field ${\mathbf{E}}_{\mathrm{DM},k}^p$ on the right-hand side instead of $-{\widehat{\mathbf{E}}}_{\mathrm{ab},k}$, and has the solution

Eq. (16)

$${\mathbf{u}}_k^p={(\mathrm{Re}\{{\mathbf{G}}_{1,k}^{†}{\mathbf{G}}_{1,k}\}+{\alpha }_{\text{probe}}^2\mathbb{I})}^{-1}\mathrm{Re}\{{\mathbf{G}}_{1,k}^{†}{\mathbf{E}}_{\mathrm{DM},k}^p\}.$$

For the experiments with SM and EFC, we applied Eq. (16) to compute probe commands for pairwise probe estimation. For experiments with AD-PSM and AD-EFC, we iteratively solved Eq. (15) using the AD-EFC framework. In all cases, we set the Tikhonov regularization parameter for probe generation as ${\alpha }_{\text{probe}}^2=0.7$; for the iterative probes, we set the optimizer tolerance to $\epsilon ={10}^{-5}$. Figure 8 shows the probe commands ${\mathbf{u}}_p$ along with the corresponding magnitude and phase of ${\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_p)$.

Fig. 8

The probe commands ${\mathbf{u}}_p$ for pairwise estimation were optimized so that the resultant dark-zone electric field ${\mathbf{E}}_{\mathrm{DM},k}({\mathbf{u}}_p)=\mathrm{sgn}\{{\mathit{\rho }}_x\}\exp \{i\pi (p-1)\mathrm{sgn}\{{\mathit{\rho }}_x\}/4\}$, as described in Sec. 3.3. The inner and outer edges of the dark zone are also shown for reference. The probe commands are close to the inverse Fourier transform of the dark zone geometry (an annulus), modulated by a horizontal sinusoid whose phase angle is proportional to the desired piston phase, and projected onto the DM actuator coordinates. Only the pupil-plane DM is modulated.

4.

Experimental Results and Discussion

We conducted a series of experiments to compare the contrast performance of AD-PSM and AD-EFC relative to SM and EFC, respectively. All experiments used an annular control region extending from $5.8{\lambda }_0/D_{LS}$ to $9.8{\lambda }_0/D_{\mathrm{LS}}$, where ${\lambda }_0=638\mathrm{nm}$ is the central wavelength of the Thorlabs MCLS1 laser diode and $D_{\mathrm{LS}}$ is the Lyot stop diameter. Each experiment consisted of 80 WFS&C iterations. For each experiment, we computed the median and 10th and 90th percentile values of the spatially averaged dark-zone contrast values for the final 50 iterations, which captured the steady-state performance of the WFS&C loop with each algorithm after converging to deepest-possible contrast. Figures 9 and 10 show example contrast vs. iteration time series with AD-PSM/SM and AD-EFC/EFC, respectively, as well as the PSF from the iteration with deepest contrast using AD-PSM and AD-EFC. In both cases, the convergence properties of AD-PSM and AD-EFC are nearly identical to their Jacobian-based counterparts, validating our proposed approach.

Fig. 9

(a) Spatially averaged dark-zone contrast versus WFS&C iteration using AD-PSM with ${\eta }_{00}=10$ and $\epsilon ={10}^{-4}$, compared to an experiment with SM. The median, 10th, and 90th percentile of the final 50 iterations, representing the steady-state behavior after convergence, are shown in green. (b) The on-axis PSF corresponding to the iteration with deepest contrast using AD-PSM.

Fig. 10

(a) Spatially averaged dark-zone contrast versus WFS&C iteration using AD-EFC with ${\alpha }_k^2={10}^{-2}$ and $\epsilon ={10}^{-4}$, compared to an experiment with EFC using the same value for ${\alpha }_k$. The median, 10th, and 90th percentile of the converged datapoints are shown in green. (b) The on-axis PSF corresponding to the iteration with deepest contrast using AD-EFC.

For SM and AD-PSM, we chose $I_{\text{target},k}=0.5{\Vert {\widehat{\mathbf{E}}}_{\mathrm{ab},k}\Vert }^2$ as the contrast target, i.e., a factor of two improvement iteration-over-iteration in spatially integrated dark zone contrast. For the Lagrange multiplier line search in SM, we let ${\mu }_k^{n+1}=1.3{\mu }_k^n$. For AD-PSM, we tested combinations of the penalty parameter ${\eta }_{00}\in \{\mathrm{10,100,1000}\}$ and the nonlinear optimization convergence tolerance $\epsilon \in \{{10}^{-2},{10}^{-3},{10}^{-4}\}$. We determined these values based on simulations and prior WFS&C experiments on HiCAT.

For EFC and AD-EFC, we selected the Tikhonov regularization matrix as ${\mathbf{\Gamma }}_k={\alpha }_k\mathbb{I}$, with ${\alpha }_k^2\in \{{10}^{-2},{10}^{-3},{10}^{-4}\}$ for the first 30 WFS&C iterations and ${\alpha }_k^2={10}^{-1}$ thereafter. We determined through previous experiments with EFC that maintaining an aggressive ${\alpha }_k$ value throughout the experiment was detrimental to the stability of the control loop after reaching the steady-state regime. For AD-EFC, we used the same set of $\epsilon$ values as AD-PSM.

For AD-EFC and AD-PSM, we compared the value of the cost function for two different starting guesses for the DM correction: ${\mathbf{a}}_k^0={\mathbf{a}}_{k-1}^{*}$, i.e., the solution of the previous WFS&C iteration, and ${\mathbf{a}}_k^0=0$. The starting guess with the lower of the two cost function values was then selected. In all cases ${\mathbf{a}}_k^0=0$ was ultimately chosen as the starting guess.

Figures 11 and 12 show the statistics of the post-convergence spatially averaged dark-zone contrast achieved with AD-PSM and AD-EFC for each combination of regularization (${\eta }_{00}$ or ${\alpha }_k$) and optimization tolerance $\epsilon$, respectively, compared to reference experiments with SM and EFC. For all parameter combinations, AD-PSM and AD-EFC equaled the contrast performance of SM and EFC, respectively. Moreover, we observed no strongly identifiable trends in achievable contrast as a function of the algorithm parameters, indicating that in all cases, AD-PSM and AD-EFC reached the contrast floor imposed by environmental instabilities. For full contrast versus iteration curves for all experiments, we refer to Fig. 13 in Appendix A.

Fig. 11

Median, 10th percentile, and 90th percentile spatially averaged contrast values achieved by AD-PSM (dark blue) as a function of optimizer tolerance, for three different values of the penalty parameter ${\eta }_{00}$, compared to SM (light blue). The contrast vs. iteration time series for the rightmost result in the left pane (${\eta }_{00}=10$, $\epsilon ={10}^{-4}$) is illustrated in Fig. 9. In all cases, the contrast performance of AD-PSM was equivalent to that of SM within statistical uncertainty.

Fig. 12

Median, 10th percentile, and 90th percentile spatially averaged contrast values achieved by AD-EFC (dark orange) as a function of optimizer tolerance, for three different values of the Tikhonov regularization parameter ${\alpha }_k$. For each ${\alpha }_k$, we also ran a reference experiment with EFC (light orange); the EFC experiment with ${\alpha }_k^2={10}^{-4}$ diverged and is not shown. The contrast versus iteration time series for the rightmost result in the left pane (${\alpha }_k^2={10}^{-2}$, $\epsilon ={10}^{-4}$) is illustrated in Fig. 10. The performance of AD-EFC was statistically equivalent to EFC for ${\alpha }_k^2={10}^{-2}$ and ${\alpha }_k^2={10}^{-3}$. As discussed in Sec. 4.1, the optimizer tolerance prevented AD-EFC from diverging in this case.

5.

Conclusions and Future Work

In this paper, we reported the first experimental demonstrations of two AD-based wavefront control algorithms, AD-PSM and AD-EFC, using the HiCAT testbed. To within statistical uncertainty, AD-PSM and AD-EFC equaled their Jacobian-based counterparts in dark-zone contrast for all combinations of parameters that we tested. These demonstrations pave the way for future experimental validation at higher contrast.

The analysis in our earlier work indicated that the largest computational gains are realized for DMs with more than $64\times 64$ actuators. The DMs currently in use on HiCAT have 34 actuators across the diameter of the active region or 952 actuators per DM in total, which is comparatively low. Therefore, our goal was to validate the fundamental capability of AD-PSM and AD-EFC to reach deep contrast, rather than to demonstrate improved computational efficiency.

In this work and in our earlier work, we utilized the L-BFGS optimization algorithm to minimize the wavefront control cost function because of its desirable convergence properties compared to first-order optimization methods and low storage requirements. Despite this, L-BFGS is known to converge slowly for poorly conditioned problems compared to methods that utilize the exact Hessian matrix.⁹ However, there exist alternative methods, such as truncated Newton algorithms, with excellent convergence properties for quadratic or nearly-quadratic cost functions, that require only the ability to evaluate Hessian-vector products, rather than the Hessian matrix itself.⁹ Hessian-vector products can be evaluated using AD in a similar fashion to gradients. Future work will explore such methods as a potentially faster approach than L-BFGS.

6.

Appendix A: Convergence Data for All Experiments

In Sec. 4, we showed the time series of spatially averaged contrast vs. WFS&C iteration for an AD-PSM with ${\eta }_{00}=10$ and $\epsilon ={10}^{-4}$ compared to an experiment with SM (Fig. 9), as well as an AD-EFC experiment with ${\alpha }_k^2={10}^{-2}$ and $\epsilon ={10}^{-4}$ compared to an EFC experiment with the same value of ${\alpha }_k^2$ (Fig. 10). We then summarized the steady-state time series statistics (median, 10th percentile, and 90th percentile) for a series of nine AD-PSM and nine AD-EFC runs with different combinations of $({\eta }_{00},\epsilon )$ and $({\alpha }_k^2,\epsilon )$, respectively (Figs. 11 and 12). In this section, Fig. 13 shows the full time series of contrast versus iteration for all 18 AD-PSM and AD-EFC experiments compared to their respective SM and EFC reference experiments.

Fig. 13

Spatially averaged dark-zone contrast versus iteration for all experiments discussed in Sec. 4 and summarized in Figs. 11 and 12. The 10th percentile, median, and 90th percentile of the steady-state contrast values (gray region) for AD-PSM and AD-EFC are shown in green. The AD-PSM experiment with $({\eta }_{00}=10,\epsilon ={10}^{-4})$ and the AD-EFC experiment with $({\alpha }_k^2={10}^{-2},\epsilon ={10}^{-4})$ are also shown in Figs. 9 and 10, respectively.

7.

Appendix B: Equivalence of Jacobian-Based Solutions and Newton’s Method

Finding solutions for EFC and SM using the Jacobian matrix is equivalent to minimizing their respective cost functions with respect to ${\mathbf{a}}_k$ using Newton’s method. Newton’s method is a second-order optimization technique that utilizes second-derivative information about the cost function, given by the local Hessian matrix ${\mathbf{H}}_k({\mathbf{a}}_k^n)$ at any point ${\mathbf{a}}_k^n$ in the DM command parameter space. Given an initial guess for the solution ${\mathbf{a}}_k^0$, the full Newton update is given as⁹

For cost functions that are exactly quadratic, including EFC and SM, Newton’s method converges in a single iteration.

For general numerical optimization problems, Newton’s method is rarely used in practice because forming the Hessian matrix explicitly is expensive. On the other hand, quasi-Newton methods, such as the BFGS algorithm or the limited-memory BFGS (L-BFGS) variant, can approximate ${\mathbf{H}}_k^{-1}$ using changes in $\partial J_k/\partial {\mathbf{a}}_k^T$ over successive optimization iterations. As a consequence, they are substantially less computationally expensive. Although quasi-Newton algorithms do not converge as rapidly as Newton’s method, and in particular can converge slowly for poorly conditioned problems, they are nonetheless superior to purely first-order methods such as steepest descent.⁹

As we show below, the Hessian matrix for EFC and SM has an analytic expression in terms of the Jacobian ${\mathbf{G}}_k$ given by ${\mathbf{H}}_k=2(\mathrm{Re}\{{\mathbf{G}}_k^{†}{\mathbf{G}}_k\}+\mathbf{C})$, where $\mathbf{C}$ is a symmetric, positive-definite matrix given by ${\mathbf{\Gamma }}_k^T{\mathbf{\Gamma }}_k$ for EFC and $\mathbb{I}/{\mu }_k^n$ for SM. Our approach is to replace this full Newton iteration, requiring computation of the Jacobian, by a series of cheaper quasi-Newton iterations instead, requiring only computation of the gradient $\partial J_k/\partial {\mathbf{a}}_k^T$, which we achieve using AD.

The cost function for the EFC algorithm described in Sec. 2 can be written in the form

Expanding Eq. (18) and recalling that ${\mathbf{E}}_{\mathrm{DZ},k}={\mathbf{G}}_k{\mathbf{a}}_k+{\widehat{\mathbf{E}}}_{\mathrm{ab},k}$, the cost function has the form

This is a linear system of equations that we can solve for the optimal correction ${\mathbf{a}}_k^{*}$:

The Hessian matrix is given as

Since both terms in ${\mathbf{H}}_k$ are positive definite, ${\mathbf{H}}_k$ is positive definite as well, confirming that the solution is a minimum of the cost function.

We will now show that the solution ${\mathbf{a}}_k^{*}$ obtained above is the same as the solution obtained by applying a single iteration of Newton’s method to the EFC cost function. Let ${\mathbf{a}}_k^0$ be an initial guess for the solution. Newton’s method produces an iterate of the form⁹

Combining Eqs. (20) and (22)

Inserting back into Eq. (23)

Inserting the definition of the Hessian matrix from Eq. (22), we see that the Newton iterate ${\mathbf{a}}_k^1$ is identical to the analytical solution in Eq. (21)

As we described earlier, the same result holds if minimizing ${\mathcal{L}}_{\mathrm{SM},k}$ in Sec. 2.1 with respect to ${\mathbf{a}}_k$.

8.

Appendix C: Fast Convolutional DM Model

Consider a DM with $N_A$ actuators along each side (i.e., $N_{\mathrm{act}}=N_A^2$) whose surface $s(x,y)$ can be modeled as a linear superposition of identical influence functions $f(x,y)$

For fixed actuator spacing along the horizontal and vertical directions, we can rewrite the above summation as a convolution between a weighted Dirac comb function and the influence function

Fourier transforming both sides transforms the convolution operation into a multiplication:

We define $\tilde{s}(f_x,f_y)\triangleq \mathcal{F}\{s(x,y)\}$ and $\tilde{f}(f_x,f_y)\triangleq \mathcal{F}\{f(x,y)\}$:

We next define the discretized surface and influence function arrays $\tilde{s}$ and $\tilde{f}$ such that

We can write this more succinctly as the following sequence of element-wise and matrix products

The term in the parentheses is more commonly referred to as the matrix Fourier transform²¹ or matrix triple product Fourier transform²⁵ of $\mathbf{A}$, which we denote as follows:

The final step is to compute an inverse discrete Fourier transform to obtain the desired discrete DM surface $\mathbf{s}$, which is carried out most efficiently using the inverse fast Fourier transform, yielding the final result

For DMs whose active actuators are a subset of the $N_A\times N_A$ grid modeled above, only the elements of $\mathbf{A}$ corresponding to active actuators are set to nonzero values.

9.

Appendix D: Adjoint Model for EFC Cost Function

In Sec. 2.2, we describe the cost function for the EFC algorithm for a single correction wavelength. Here, we derive its RMAD adjoint model, which computes the derivative $\partial J_{\mathrm{EFC},k}/\partial {\mathbf{E}}_{\mathrm{DM},k}$.

We begin by writing the cost function as a series of operations evaluated sequentially:

We now apply the RMAD gradient rules^7,26 to each step in reverse order to derive the adjoint model, letting $\overline{\mathbf{x}}\triangleq \partial J_k/\partial {\mathbf{x}}^T$ for any variable $\mathbf{x}$:

Combining and simplifying, we see that the desired gradient is given as

This gradient is passed to the next block of the adjoint model, which in this case is the adjoint model for propagation through the coronagraph, which, referring to Fig. 4, evaluates the derivatives of $J_{\mathrm{EFC},k}$ with respect to the surfaces ${\mathbf{s}}_{1,k}$ and ${\mathbf{s}}_{2,k}$, respectively. We refer to our earlier work for a derivation of the coronagraph propagation adjoint model.⁷