Iterative X-ray spectroscopic ptychography1

Spectroscopic ptychography is a powerful technique to determine the chemical composition of a sample with high spatial resolution. This paper presents a novel algorithm to iteratively solve the spectroscopic blind ptychography problem.


Introduction
X-ray spectro-microscopy is a powerful technique to study the chemical and morphological structure of a material at high resolution. The contrast of the material under study is recorded as a function of photon energy, and this spectral absorption contrast can later be used to reveal details about its chemical, orbital or magnetic state (Stö hr, 2013;Koningsberger & Prins, 1988). The idea is that, because different chemical components interact differently with the beam at different energies, the composition map of a sample can be solved by using measured reference spectra (a dictionary).
Compared with standard lens-based microscopy, X-ray ptychography can provide much finer spatial resolution, while also providing additional phase contrast of the sample (Nellist et al., 1995;Chapman, 1996;Rodenburg & Faulkner, 2004;Rodenburg et al., 2007). Ptychography is based on retrieving the phase of diffraction data recorded to a numerical aperture that is far larger than what X-ray optics can technically achieve. In ptychography, the probe (illumination) is almost never completely known, so a joint recovery problem (sample and probe) is typically considered, referred to as blind ptychography. Several algorithms to solve both standard and blind ptychography problems have been published in the literature, which also consider a variety of additional ISSN 1600-5767 experimental challenges (Maiden & Rodenburg, 2009;Thibault et al., 2009;Thibault & Guizar-Sicairos, 2012;Wen et al., 2012;Marchesini et al., 2013;Horstmeyer et al., 2015;Hesse et al., 2015;Odstrci et al., 2018;Chang et al., 2019a).
As in standard spectro-microscopy, it is possible to perform spectroscopic ptychography by recording diffraction data at different X-ray photon energies. In recent years, spectroptychography has become an increasingly popular chemical analysis technique (Beckers et al., 2011;Maiden et al., 2013;Hoppe et al., 2013;Shapiro et al., 2014;Farmand et al., 2017;Shi et al., 2016). However, the standard methodology involves independent ptychographic reconstructions for each energy, followed by component analysis, i.e. spectral imaging analysis based on a known reference spectrum or multivariate analysis (Adams et al., 1986;Lerotic et al., 2004;Shapiro et al., 2014;Yu et al., 2018). More recently, a low-rank constraint (Vaswani et al., 2017) for multi-channel samples was proposed, together with a gradient descent algorithm with spectral initialization to recover the higher-dimension phase-retrieval problem (without component analysis). Other work proposed a hierarchical model with a Gaussian-Wishart hierarchical prior and developed a variational expectation-maximization algorithm (Liu et al., 2019). Also, a matrix-decomposition-based lowrank prior (Chen et al., 2018) has been exploited to reconstruct dynamic time-varying targets in Fourier ptychographic imaging.
In this paper we propose a novel technique to solve the blind X-ray spectro-ptychography problem, based on coupling the diffraction data from each photon energy and iteratively retrieving the chemical map of the sample. The proposed algorithm, referred to as SPA (spectroscopic ptychography with ADMM), works with both completely and partially known reference spectra. The method is designed using the alternating direction method of multipliers (ADMM) (Glowinski & Tallec, 1989;Chang et al., 2019a) framework, employing also total variation (TV) regularization (Rudin et al., 1992) on the chemical map. Compared with the standard two-step methods, the proposed joint reconstruction algorithm can generate much higher quality results without presenting the phase ambiguity problem inherent to two-step methods. 21 The simulation analysis shows the efficient convergence ratio of SPA and demonstrates the increased robustness of the method to large step sizes, being able to retrieve features lost when using standard two-step methods. The algorithm is described and analyzed with and without TV regularization for both partially and completely known dictionary cases.

Spectroscopic ptychography model
The main operators used in this section are given in Table 1.
Given L different energies of X-rays going through a sample illuminated by a probe ! 2 C m , a collection of phaseless intensities fI l g LÀ1 l¼0 are measured in the far field, such that with Poisson fluctuation caused by photon counting we have Here Y ¼ ðY 0 ; Y 1 ; . . . ; Y LÀ1 Þ 2 C N;L is the sample contrast map for each X-ray energy, Aðw; ÁÞ is the forward operator for ptychography for a given probe w, Poi denotes the Poissonnoise contamination, and the notations j Á j; ðÁÞ 2 denote the pointwise absolute and square values of the vector, respectively. Note that the probe ! and each column of contrast maps fY l g are all 2D images, written as vectors by a lexicographical order. The relationship between the contrast map Y l observed by ptychography at each energy and an unknown sample elemental map X made of C elements is governed by the spectral contrast of each element, stored in a 'dictionary' D of known values. Specifically, following similar notation to Chang, Enfedaque et al. (2018), the bilinear operator A : C m Â C N ! C m is defined as Að!; uÞ :¼ n ½F ð! S 0 uÞ T ; ½F ð! S 1 uÞ T ; . . . ; where F denotes the discrete Fourier transform, denotes the Hadamard product (pointwise multiplication) of two vectors, and S j 2 R m;m is a binary matrix that defines a small window with the index j and size m over the entire image u (taking small patches out of the entire image). For different energies, assuming that a spectrum dictionary D 2 C C;L (or its absorption part) is measured in advance, having C components for different materials or particles, and given a sufficiently thin specimen, the sample contrast maps can be approximated by first-order Taylor expansion expðXDÞ ' 1 þ XD 31 as with X ¼ ðX 0 ; X 1 ; . . . ; X CÀ1 Þ 2 R N;C þ being the elemental thickness map of the sample (each column of the thickness map denotes the thickness of each component in the object).  Table 1 Main variables and operators defined in Section 2.

Notations
Explanations Measured intensities X 2 R N;C þ Elemental thickness maps of the sample D 2 C C;L Spectrum dictionary Y 2 C N;L Sample spectral contrast maps ! 2 C m Probe S j 2 R m;m Binary matrix to take image patches Forward operator for ptychography G½Aðw; uÞ; I l Poisson likelihood estimation D r 2 R C;L Real part of spectrum dictionary X Non-negative thickness constraint W Normalized constraint of the probe TV Total variation regularization To determine the thickness map X, with a completely known spectrum D, one has to solve the following problem: To find X and !; such that with non-negative thickness constraint set X ¼ fX ¼ ðX n;c Þ 2 R N;C : X n;c ! 0; 0 n N À 1; 0 c C À 1g. Letting the illumination be normalized, i.e. ! 2 W :¼ f! 2 C m : k!k ¼ 1g, the total variation regularized nonlinear optimization model can be established by assuming the piecewise smoothness of the thickness map as SP : min . . . ; f NÀ1 Þ T 2 R N , derived from the maximum likelihood estimate of Poisson-noised data , TV denotes the standard total variation seminorm (Rudin et al., 1992) to enforce the piecewise smooth structure of X c [the cth column of the mixing matrix (thickness map) defined in equation (3)], and is a positive constant to balance the regularization and fitting terms (larger produces stronger smoothness). I X ðXÞ, We remark that this is a convenient way to enforce hard constraints within an optimization formulation.
Experimentally, only the real-valued part (absorption) of the dictionary D r is measured. As X is real valued, we consider the following relation: where D r :¼ <ðDÞ. Similarly, we derive the following for spectroscopic ptychography with an incomplete dictionary (SPi): Remark: Rather than solving the ptychography imaging independently for each energy, we use the low-rank structure of the recovery results of different energies, i.e. the rank of the matrix Y À 1 is no greater than that of X.

Proposed iterative algorithm
ADMM (Glowinski & Tallec, 1989) is a powerful and flexible tool that has already been applied to both ptychography (Wen et al., 2012;Chang et al., 2019a) and phase tomography problems (Chang et al., 2019b;Aslan et al., 2019). In this work we also adopt the ADMM framework to design an iterative joint spectro-ptychography solution. We construct the proposed algorithm considering both complete and incomplete dictionary cases.

Complete dictionary
On the basis of the spectro-ptychography model [equation (5)] for a complete dictionary of spectra, we design the proposed SPA algorithm as described below.
Let DD Ã 2 C C;C be non-singular, where D Ã denotes the Hermitian matrix of D, i.e. D Ã :¼ conjðD T Þ. Considering the constraint in equation (3), the following equivalent form can be derived: withD D :¼ D Ã ðDD Ã Þ À1 2 C L;C . Accordingly, the following equivalent model can be considered, by introducing auxiliary variables fZ l g: The benefits of considering equation (8) instead of equation (3) lie in the fact that (i) the multiplier will be a low-dimensional variable, since the dimension of Y is much higher than that of X, and (ii) the subproblem with respect to the variable X can be more easily solved. An equivalent saddle point problem for equation (9), based on the augmented Lagrangian, can be derived as max Ã;À min !;X;Y;Z L ; ð!; X; Y; Z; Ã; ÀÞ with the multipliers Ã :¼ ðÃ 0 ; . . . ; Ã LÀ1 Þ and À, where hÁ; Ái denotes the inner product of two vectors (or trace norms for two matrices). The above saddle point problem can be solved by alternating minimization and update of the multipliers. We first define each sub-minimization problem. The ! subproblem, with the additional proximal term, can be expressed as The first-order gradient of the above least-squares problem (without constraint) is given as Consequently, the projected gradient descent scheme with preconditioning can be derived as research papers with parameter 1 > 0 in order to avoid division by zeros, and Proj W ð!Þ :¼ !=k!k. Here the parameter 1 is heuristically set to be a small scalar related to the maximum value of P l;j jS j Y l j 2 , e.g. 1 ¼ 0:1k P l;j jS j Y k l j 2 k 1 . The X subproblem can be expressed as where ðÁÞ c denotes the cth column of a matrix. Since it is common practice to solve the total variation denoising problem by using a first-order operator-splitting algorithm (Wu et al., 2011;Chambolle & Pock, 2011), we directly give the approximate solution below: with Denoise ðu 0 Þ :¼ arg min u TVðuÞ þ 1 2 ku À u 0 k 2 . Here we remark that, to seek the exact solution with this positivity constraint, one may need more auxiliary variables and inner loops (Chan et al., 2013). For simplicity, we did not exactly solve the constraint problem, and instead, the above approximation is derived by the standard TV-L2 denoising without constraint and then a projection to the positivity constraint set.
The Y subproblem, with additional proximal term ð 2 =2ÞkY À Y 0 k 2 and previous iterative solution Y 0 , is expressed as where 2 is a positive scalar similar to the parameter 1 . By calculating the first-order gradient of the above leastsquares problem, one has with identity operator I, where Q :¼ ðQ 0 ; Q 1 ; . . . ; Q LÀ1 Þ 2 C N;L with Q l :¼ P j S T j ½! Ã F Ã ðÃ l þ Z l Þ. Equation (17) is actually the Sylvester equation (Sylvester, 1884;Simoncini, 2016).
Assuming that the positive HermitianD DD D Ã has the singular value decomposition (SVD)D DD D Ã ¼ VS V Ã , with diagonal matrix (diagonal elements are singular values) S 2 R L;L and unitary matrix V 2 C L;L , and by introducingŶ Y : such that the closed-form solution can be expressed as For the Z subproblem, we have ) which gives withẐ Z l :¼ Að!; Y l Þ À Ã l . The above calculations and the update of the multipliers form the basis of the baseline SPA algorithm, summarized in Appendix A.

Incomplete dictionary
A complete dictionary is often difficult to obtain without an independent experiment prior to a spectro-ptychography experiment. The material's components and their chemical states are often not known in advance. Moreover, the real part of the refractive index component is often not well known (Henke et al., 1993). It is more difficult to measure because it requires interferometric or reflectometry measurements rather than simple absorption spectroscopy measurements, and reflectometry experiments are less commonly done. While the Kramers-Kronig relationships relate real and imaginary parts, the relationship requires a spectral measurement from 0 to infinity, which is not possible to measure in finite time. Standard techniques to extend absorption spectra can only produce approximate values in the imaginary component. Hence, it is attractive in practice to provide a version working with the real part only.
In this subsection, we propose a variation of the SPA algorithm to solve the joint spectro-ptychography problem when the dictionary of spectra is only partially know, based on the model proposed in equation (7). By assuming that D r has full row rank, i.e. D r D Ã r is non-singular, withD D r :¼ D T r ðD r D T r Þ À1 known in advance, we have research papers Consequently, the following equivalent problem can be solved instead of equation (7): Similarly to the previous subsection, introducing the multiplier À r and auxiliary variable Z yields the saddle point problem below, with the help of the augmented Lagrangian of equation (24): Below, we focus only on the differences with respect to Algorithm 1 (Appendix A). For the X subproblem, we have X ? :¼ arg minð=Þ P c TVðX c Þ þ I X ðXÞ þ 1 2 kX À ½<ðY À 1ÞD D r À À r k 2 : Hence we get For the Y subproblem with proximal terms kY À Y 0 k 2 , we have which results in the following equations with respect to the real and imaginary parts, respectively: Then, the real part of Y can be solved by equations (19) and (20), while the imaginary part can be simply computed by The overall SPA algorithm with an incomplete dictionary is summarized in Appendix B.

Figure 2
Spectrum dictionaries [three different materials m0 (PMMA), m1 (PS) and m2 (constant)], with real part (a) and imaginary part (b). The x and y axes denote the order of the ten spectra and different energies, respectively.

Simulation and reconstruction results
In the simulation analysis of the proposed algorithms we consider the synthetic thickness maps of three different materials, extracted from three (RGB) channels of a natural color image (after thresholding and shift, consisting of 256 Â 256 pixels), shown in Fig. 1. The real part of the spectrum dictionary [for two different materials, (a) PMMA (polymethyl methacrylate) and (b) PS (polystyrene), plus (c) a constant with respect to ten different energies] was measured at the Advanced Light Source (Yan et al., 2013), and the imaginary part was derived using the Kramers-Kronig relations (Kronig, 1926). Both real and imaginary part dictionaries are shown in Fig. 2.
The ptychography measurements are simulated with Poisson noise contamination, using a single grid scan at each energy. A standard zone plate with annular shape diffracts an illumination (Fig. 3a) onto the sample after the focused probe (Fig. 3b) has gone through an order-sorting aperture. The zone plate annular aperture is mapped onto the detector by geometric magnification as 'outer diameter' (in mm) and corresponds to an annular ring on the detector of dimension (outer diameter /detector pixel size) Â (detector distance / focal distance). The illumination probe [ Fig. 3(b)] has a beam width (FWHM) of 16 pixels. The relationship between pixels and actual dimensions in the far-field approximation is as follows: illumination pixel (real-space) dimensions = (wavelength Â detector distance) /(detector number of pixels Â detector pixel size). The zone plate's distance from the sample is assumed to be adjusted proportionally with energy to keep the sample in focus, as is usually done experimentally. We also

Figure 4
Reconstruction results using a known dictionary of spectra from Poisson-noised data with SNR = 29.2 dB and scan step size = 32. (a)-(c) Standard twostep method; (d)-( f ) SPA without regularization; (g)-(i) SPA with TV regularization. The recovered probes are shown in the right column for the twostep method, SPA and SPA with TV (from top to bottom). assume that the detector distance is adjusted to maintain the spatial frequencies on the same detector pixels.
In order to evaluate the recovered results, the signal-tonoise ratio (SNR) in dB is used, which is defined below: SNRðX; X g Þ ¼ À10 log 10 kX À X g k 2 =kXk 2 ; where X g corresponds to the ground-truth thickness. We compare the proposed iterative SPA algorithm with the standard two-step method. The two-step method consists of (i) performing ptychography reconstruction using a joint illumination, then (ii) performing spectroscopy analysis with a known dictionary (or known real part), and finally (iii) correcting the phase ambiguity for different energies.
When assessing the performance of SPA, we consider both with and without regularization cases, where we simply set the regularization parameter ¼ 0 and slightly adjust the algorithm by replacing Step 2 with for baseline SPA and X kþ1 ¼ maxf0; <ðY k À 1ÞD D r À À k r g ð 33Þ for the incomplete dictionary case.

Reconstruction quality
The first simulation assesses the reconstruction quality achieved by the proposed SPA algorithm, compared with the two-step method, when using a scan step size of 32 pixels. Figs. 4 and 5 depict the reconstructed images when using complete and incomplete dictionaries, respectively. The SPA simulations are performed without and with regularization in rows 2 and 3, respectively, of Figs. 4 and 5. Visually, we can see obvious artifacts in the recovered images when using the twostep method (first row of Figs. 4 and 5). Such artifacts are greatly enhanced when the image is reconstructed using SPA. Specifically, clear improvements can be identified in the regions corresponding to the red and blue circles for all three materials in both Fig. 4 and Fig. 5. The SNRs of the recovery results parallel the visual analysis. For the completely known dictionary, the two-step method achieves an SNR of 14.0 dB for the above simulation, whereas SPA achieves 18.1 dB (no regularization) and 18.8 dB (regularization). In the partially Reconstruction results using a partially known dictionary of spectra from Poisson-noised data with SNR = 29.2 dB and scan step size = 32. (a)-(c) Standard two-step method; (d)-( f ) SPA without regularization; (g)-(i) SPA with TV regularization. The recovered probes are shown in the right column for the two-step method, SPA and SPA with TV (from top to bottom). known dictionary simulation, the SNRs are 13.8, 15.8 and 16.7 dB, for the two-step method and SPA with no regularization and regularization, respectively, which achieves a comparative gain of more than 2 dB, similarly to the known dictionary case.
The phase ambiguity is an inherent problem of the two-step method that causes a loss in reconstruction accuracy. For example, for the simulation shown in Fig. 4, the SNR without phase correction is only 12.3 dB, reaching 14.0 dB after applying correction. Even when using an effective phase correction post-process, SPA proves to be more efficient for the simulations performed: higher-quality reconstructions are achieved overall, and there is no need to correct the phase ambiguity because of the iterative reconstruction exploiting the low-rank structure and positivity constraint of the thickness function.

Robustness and convergence
The following simulation assesses the robustness of the proposed algorithm when varying the scanning step sizes. The quantitative results of this simulation are presented in Table 2.
The results demonstrate the enhanced robustness of SPA when handling larger step sizes, compared with the reference two-step method, achieving up to 10 dB increase in SNR. To permit a better visual analysis, we provide the reconstruction results of the three algorithms with 40 pixels step size in Fig. 6. The figure highlights the dramatic improvement achieved by SPA compared with the standard two-step method when reconstructing low-redundancy ptychographic data. Specifically, we can see how the features within the blue and red circles are almost lost in the two-step reconstruction, while they can be clearly observed when reconstructing using SPA.
Generally speaking, to make the proposed algorithms work, the basic condition is to assume D has full row rank such that Reconstruction results using a known dictionary of spectra from Poisson-noised data with SNR = 29.0 dB and scan step size = 40. (a)-(c) Standard twostep method; (d)-( f ) SPA without regularization; (g)-(i) SPA with TV regularization. The recovered probes are shown in the right column for the twostep method, SPA and SPA with TV (from top to bottom). DD Ã is non-singular. However, the performance should also rely on the similarity of spectral elements. Here we introduce a factor s 2 ½0; 0:5 to generate a new dictionary D s 2 C 3Â10 with D s ½1; k ¼ ð1 À sÞD½1; k þ sD½2; k, D s ½2; k ¼ ð1 À sÞD½2; k þ sD½1; k 8 1 k 10. We know that (1) D 0 ¼ D and (2) the first two rows are exactly the same if s ¼ 1=2 (D 1=2 does not have full row rank). (See Fig. 7 for the dictionaries with s = 0.1, 0.3 and 0.45.) Therefore, the parameter s can be used to control the similarity of the new spectral dictionary (larger s implies higher similarity). We test the impact of the proposed SPA algorithm by the different similarity of spectral dictionaries (see the SNR changes in Fig. 8  Synthetic spectrum dictionaries D s for different s [three different materials m0 (PMMA), m1 (PS) and m2 (constant)], for real part (left) and imaginary part (right). The x and y axes denote the spectrum and different energies, respectively.

Figure 8
SNR changes versus similarity parameter s.
the spectral dictionaries become similar (the parameter s gets close to 0.5). Hence, to get a better reconstruction, one should design the experiments with little similarity in spectral dictionaries. The last simulation depicts the error curve achieved by the SPA algorithm, shown in Fig. 10. The results demonstrate a steady decrease of the successive errors of the proposed algorithm, both with and without regularization.

Conclusions
This paper presents the first iterative spectroscopy ptychography solution. The proposed SPA algorithm is based on a novel spectro-ptychography model and it is constructed considering both a completely known and a partially known dictionary. Numerical simulations show that SPA produces more accurate results with clearer features compared with the standard two-step method. In the future, we will extend our work to thicker samples, where the first-order Taylor expansion is not sufficiently accurate. We also plan to investigate the use of Kramers-Kronig relationships (Hirose et al., 2017), explore the case using a completely unknown dictionary and further provide software for real experimental data analysis.
APPENDIX A Algorithm 1: SPA  Reconstruction results by proposed SPA using different dictionaries of spectra (s = 0.1, 0.3, 0.45), from Poisson-noised data with SNR = 29.2 dB and scan step size = 32.

Figure 10
Error kx k À X kÀ1 k=kX k k variation versus iteration number for SPA without (a) and with TV regularization (b). The x and y axes denote iteration numbers and errors, respectively.
(1) Update the probe ! kþ1 by the one-step projected gradient descent method: 1 ! k ð P l;j jS j Y k l j 2 Þ þ 1 1 þ P l;j F Ã ðZ k l;j þ Ã k l;j Þ S j ðY k l Þ Ã ð P l;j jS j Y k l j 2 Þ þ 1 1 # ; with 1 ¼ 0:1k P l;j jS j Y k l j 2 k 1 .