# Norm-Minimized Scattering Data from Intensity Spectra

Authors: Alexander Seel; Arman Davtyan; Ullrich Pietsch; Otmar Loffeld
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We apply thel1 minimizing technique of compressive sensing (CS) to nonlinear quadratic observations. For the example of coherent X-ray scattering we provide the formulas for a Kalman filter approach to quadratic CS and show how to reconstruct the scattering data from their spatial intensity distribution.


---

## Body

## 1. Introduction

The rapidly growing Si technology in semiconductor electronics [1, 2] opens the possibility to grow III–V inorganic nanowires, such as GaAs, InAs, or InP, which were supposed to have potential [3, 4] to become building blocks in a variety of nanowire-based nanoelectronic devices, for example, in nanolaser sources [5] or nanoelectronics [6]. Such epitaxially grown nanowires are repeating the crystal orientation of the substrate and usually grow in Wurtzite (WZ) or Zinc-Blende (ZB) structure differing in the stacking sequences ABABAB and ABCABCABC, respectively, of the atomic bilayers. Theoretical predictions on the electronic properties [7] of these nanowires show that stacking sequences with WZ and ZB segments considerably differ in the conductivity. However, during the nanowire growth stacking faults, the mixing of ZB and WZ segments takes place, and twin defects [8] appear. As these defects have their own impact on the conductivity and band structure there is great interest in knowing the exact stacking sequence which can be studied by, for example, Transmission Electron Microscopy [9]. But, as this is a destructive method, it is impossible to use the nanowire after the structural studies. Nowadays the 3rd generation synchrotron sources and rapidly developing focusing devices like Fresnel Zone Plates opens new fields of nondestructive X-ray imaging. For example, in the Coherent X-ray Diffraction Imaging (CXDI) experiments, an isolated nanoobject is illuminated with coherent X-ray radiation and the scattered intensity is measured by a 2D detector [10, 11] under the Bragg condition. The diffraction patterns are structure-specific and encode the information about the electron density of the sample and thus the stacking sequence of the atomic bilayers formally by Fourier transform. However, because of sensor physics, the phase information is lost in CXDI measurement since the measured intensity pattern is the modulus square of the scattered X-rays and no inverse Fourier transform can be directly applied to recover the stacking sequence. The classical approach with the Patterson function [12] fails as the number of expected randomly distributed stacking faults is too high. Other inversion algorithms [13, 14] could be used instead to reconstruct the lost phase information: although dual space iterative algorithms [15] have been shown to converge under specific conditions [16] there still remain some convergence problems for a number of cases when not enough preliminary information regarding the structure of the object is previously known [17]. Indeed ptychography type of experiments [18, 19] were suggested to determine at least the relative phases of the bilayers by interfering adjacent Bragg reflections. But this type of experiments is more difficult to realize, as it requires higher stability in comparison to conventional CXDI and the longer measurement time can, however, influence the sample [20].Due to the principal loss in phase information the scattering data can be considered to be undersampled and algorithms tailored to this lack of information, like basis pursuit approaches [21] realized by minimizing the l1 norm of sparse vectors in an appropriate basis or phase retrieval via matrix completion [22–24], could be tested for reconstruction from conventional CXDI measurements: as these data are recorded without structured illumination for applying matrix completion we focus on utilizing the l1 minimizing technique to reconstruct undersampled sparse signals [25, 26] as vectors from linear mappings between the signal and the observations. For an overview of this so-called compressive sensing (CS) see the textbook [27]. The method of CS aims at a signal’s sparse support which can be estimated by, for example, Kalman filtering [28]. As the underlying filter formulas also relate observations and signals to be reconstructed by linear mappings this technique meets CS and was also used to explicitly minimize the l1 norm by so-called pseudomeasurements [29]. For an example see the reconstruction from a random sample of Fourier coefficients [30], where the dimensionality of the underlying Kalman gain matrices can be reduced utilizing the null space of the sensing matrix [31]. Even for nonlinear mappings of signals Kalman filtering applies by using Jacobians rather than constant sensing matrices. These extended Kalman filters (EKF) are used for, for example, tracking issues [32] or robotics [33] and matching the sensing problem of the quadratic nonlinearities in the spatial intensity distribution of CXDI.The paper is organized as follows. In Section2 we set up the observation model for modulus-squared amplitudes of intensity distributions in coherent X-ray scattering and point out the relation to the l1 minimization. In Section 3 we give a brief overview of the linear Kalman filter model and show how to incorporate the complex nonanalytic l1 norm as a linearized observation to meet the minimizing strategy of CS. In this framework we prove a convergence concept for the l1 minimization in the reconstruction scheme. In Section 4 we apply our findings to 1D and 2D Fourier data and remark in Section 5 on future considerations.


## 2. Motivation

### 2.1. Intensity Spectra in Coherent Scattering Experiments

Exposing crystals to nondestructive coherent X-ray radiation ind dimensions the amplitude of the elastically scattered radiation is proportional to the Fourier transform [34](1)Sq=∫ddxeiq∣rbrof the electron density b(r), where the vector q≔kout-kin is a parametrization of the direction where radiation is detected and (multiplied by Planck’s constant) also describes the momentum transfer in a kinematical scattering theory; kin is the incident direction whereas kout represents the outgoing direction for radiation with wave vectors kout=kin=2π/λ of wavelength λ (see Figure 1).Figure 1Following standard textbooks, for example, [34], one deals with two sets of basis vectors to characterize the scattering. The spatial vectors are represented in the basis a1,…,ad of the grid whereas the reciprocal basis b1,…,bd is used for wave vectors q,kin,kout relying on the normalization(2)aj∣bk=2πδjk,j,k=1,…,d.As the units of lengths and wave vectors are carried by the basis sets the scalar products read, with dimensionless factorsyj and κk,(3)r=∑j=1dyjaj,q=∑k=1dκkbk,q∣r=2π∑j=1dyjκj.Restricting to a finite grid withn1,n2,…,nd lattice sites in each direction with periodic boundary conditions the possible grid positions and wave vectors allowing for a discrete Fourier transform1 read, for all j=1,2,…,d,(4)yj∈0,1,…,nj-1,κj∈0nj,1nj,…,nj-1nj.Thus the scalar product separates into thed spatial dimensions according to(5)q∣r=∑j=1d2πkjrjnj,kj,rj∈0,1,…,nj-1,where kj are related to the lattice positions in direction of the lattice constant aj and rj refer to a discretized wave vector in direction of the corresponding reciprocal basis bj. For example, a d=2-dimensional regular hexagonal lattice a1,a2 encloses 120∘angles (see Figure 2).Figure 2The bold parallelogram represents periodic boundary conditions with, for example,n1=4 and n2=3 fragmenting the 1st Brillouin zone, the elementary cell spanned by b1,b2 in reciprocal space, into a subgrid according to (4). With the normalization (2) the reciprocal lattice encloses angles of 60∘ with a fixed orientation with respect to a1,a2.


### 2.2. Setting

Of course for the application the formulas are only needed up tod=3 dimensions: nanowires can be characterized by the stacking sequence of atomic bilayers which are shifted laterally and vertically with respect to each other. In coherent X-ray diffraction summing up all scattered radiation from a certain bilayer perpendicular to the growth direction a3 yields a complex scattering amplitude xk∈C associated with the kth bilayer spanned by a1,a2. Due to the hexagonal lattice and with respect to an arbitrary reference bilayer there are three different phase factors:(6)1,exp⁡2πi32κ1+κ2,exp⁡2πi3κ1+2κ2for the amplitudes left [34, 35]. Wave vectors q sensitive to the arrangement of the bilayers are selected from directions related to the Bragg condition by κ1-κ2≠3N with N,κ1,κ2∈Z yielding 1,exp(-2πi/3),exp(2πi/3), which directly relates to the three possible lateral shifts. With respect to the experimental setup [17] this is satisfied by κ1=0 and κ2=1. Recovering the stacking sequence of these relative phase factors by varying κ3 directly reveals the Zinc-Blende or Wurtzite structure of the wire along with their stacking faults.Carrying out the remaining summation over all equidistant bilayers in growth directiona3 mathematically corresponds to the 1D discrete Fourier transform(7)Sr=∑k=0n-1e2πikr/nxk,r=0,1,…,n-1 of the periodically continued complex scattering amplitudes xk combined with the vector x. Each nonzero entry represents one bilayer and r refers to the discretized component κ3 of the q vector in the reciprocal basis corresponding to the growth direction. As the outcome of the detector is the measured intensity rather than the amplitude the observation is the squared signal (8)Sr2=xTrx,r=0,1,…,n-1,where the Fourier coefficients build a Hermitian Toeplitz matrix(9)Tr≔e-2πir/np-qpq∈Cn×n,r=0,1,…,n-1,showing orthogonality TrTs=nδrsTr with respect to multiplication and a completeness property T0+T1+⋯+Tn-1=n1n with respect to summation. Clearly, due to the squaring, the signal (8) is invariant under both the transformation x→eiϕx with any global phase ϕ∈R and the reflection xk→x¯n-1-k,k=0,…,n-1, on the lattice. As the spatial directions exhibit periodic boundary conditions the signal also remains invariant under discrete translations xk→xk+q,q∈Z, on the grid.On the contrary, setting theq component κ3 corresponding to the growth direction a3 in the reciprocal basis to zero, (1) examines the 2D structure of all the M bilayers simultaneously. If they are assumed to be identical without lateral2 shifts the scattering amplitudes of the bilayer’s lattice sites are 2D real data sets Xk1k2 with the discrete Fourier transform(10)
S
r
1
r
2
=
M
∑
k
1
=
0
n
1
-
1
∑
k
2
=
0
n
2
-
1
exp
⁡
2
π
i
k
1
r
1
n
1
exp
⁡
2
π
i
k
2
r
2
n
2
X
k
1
k
2for fixed r1=0,1,…,n1-1 and r2=0,1,…,n2-1 discretizing κ1 and κ2. Squaring (10) yields Toeplitz matrices with Kronecker product structure (A.7) for the detected signal(11)Sr1r22=M2·xTr1⊗Tr2xwith multiple row and column indices (p1p2) and (q1q2), respectively, referring to the ordinary rows p1,p2 and columns q1,q2 of the Toeplitz matrices (9):(12)Tr1⊗Tr2p1p2,q1q2=Tr1p1q1Tr2p2q2,xq1q2≔Xq1q2.The scattering amplitudes in Fourier domain can be assigned with arbitrary phases. As this leaves the given intensity distributionS2 unchanged an related infinite set of broadened scattering data is formed by convolution. But if the signal S2 is assumed to be oversampled the original scattering data in this set appear to be sparse and can be recovered by an l1 minimization run on, for example, (8) or (11)—up to translations, reflections, or global phases.


## 2.1. Intensity Spectra in Coherent Scattering Experiments

Exposing crystals to nondestructive coherent X-ray radiation ind dimensions the amplitude of the elastically scattered radiation is proportional to the Fourier transform [34](1)Sq=∫ddxeiq∣rbrof the electron density b(r), where the vector q≔kout-kin is a parametrization of the direction where radiation is detected and (multiplied by Planck’s constant) also describes the momentum transfer in a kinematical scattering theory; kin is the incident direction whereas kout represents the outgoing direction for radiation with wave vectors kout=kin=2π/λ of wavelength λ (see Figure 1).Figure 1Following standard textbooks, for example, [34], one deals with two sets of basis vectors to characterize the scattering. The spatial vectors are represented in the basis a1,…,ad of the grid whereas the reciprocal basis b1,…,bd is used for wave vectors q,kin,kout relying on the normalization(2)aj∣bk=2πδjk,j,k=1,…,d.As the units of lengths and wave vectors are carried by the basis sets the scalar products read, with dimensionless factorsyj and κk,(3)r=∑j=1dyjaj,q=∑k=1dκkbk,q∣r=2π∑j=1dyjκj.Restricting to a finite grid withn1,n2,…,nd lattice sites in each direction with periodic boundary conditions the possible grid positions and wave vectors allowing for a discrete Fourier transform1 read, for all j=1,2,…,d,(4)yj∈0,1,…,nj-1,κj∈0nj,1nj,…,nj-1nj.Thus the scalar product separates into thed spatial dimensions according to(5)q∣r=∑j=1d2πkjrjnj,kj,rj∈0,1,…,nj-1,where kj are related to the lattice positions in direction of the lattice constant aj and rj refer to a discretized wave vector in direction of the corresponding reciprocal basis bj. For example, a d=2-dimensional regular hexagonal lattice a1,a2 encloses 120∘angles (see Figure 2).Figure 2The bold parallelogram represents periodic boundary conditions with, for example,n1=4 and n2=3 fragmenting the 1st Brillouin zone, the elementary cell spanned by b1,b2 in reciprocal space, into a subgrid according to (4). With the normalization (2) the reciprocal lattice encloses angles of 60∘ with a fixed orientation with respect to a1,a2.


## 2.2. Setting

Of course for the application the formulas are only needed up tod=3 dimensions: nanowires can be characterized by the stacking sequence of atomic bilayers which are shifted laterally and vertically with respect to each other. In coherent X-ray diffraction summing up all scattered radiation from a certain bilayer perpendicular to the growth direction a3 yields a complex scattering amplitude xk∈C associated with the kth bilayer spanned by a1,a2. Due to the hexagonal lattice and with respect to an arbitrary reference bilayer there are three different phase factors:(6)1,exp⁡2πi32κ1+κ2,exp⁡2πi3κ1+2κ2for the amplitudes left [34, 35]. Wave vectors q sensitive to the arrangement of the bilayers are selected from directions related to the Bragg condition by κ1-κ2≠3N with N,κ1,κ2∈Z yielding 1,exp(-2πi/3),exp(2πi/3), which directly relates to the three possible lateral shifts. With respect to the experimental setup [17] this is satisfied by κ1=0 and κ2=1. Recovering the stacking sequence of these relative phase factors by varying κ3 directly reveals the Zinc-Blende or Wurtzite structure of the wire along with their stacking faults.Carrying out the remaining summation over all equidistant bilayers in growth directiona3 mathematically corresponds to the 1D discrete Fourier transform(7)Sr=∑k=0n-1e2πikr/nxk,r=0,1,…,n-1 of the periodically continued complex scattering amplitudes xk combined with the vector x. Each nonzero entry represents one bilayer and r refers to the discretized component κ3 of the q vector in the reciprocal basis corresponding to the growth direction. As the outcome of the detector is the measured intensity rather than the amplitude the observation is the squared signal (8)Sr2=xTrx,r=0,1,…,n-1,where the Fourier coefficients build a Hermitian Toeplitz matrix(9)Tr≔e-2πir/np-qpq∈Cn×n,r=0,1,…,n-1,showing orthogonality TrTs=nδrsTr with respect to multiplication and a completeness property T0+T1+⋯+Tn-1=n1n with respect to summation. Clearly, due to the squaring, the signal (8) is invariant under both the transformation x→eiϕx with any global phase ϕ∈R and the reflection xk→x¯n-1-k,k=0,…,n-1, on the lattice. As the spatial directions exhibit periodic boundary conditions the signal also remains invariant under discrete translations xk→xk+q,q∈Z, on the grid.On the contrary, setting theq component κ3 corresponding to the growth direction a3 in the reciprocal basis to zero, (1) examines the 2D structure of all the M bilayers simultaneously. If they are assumed to be identical without lateral2 shifts the scattering amplitudes of the bilayer’s lattice sites are 2D real data sets Xk1k2 with the discrete Fourier transform(10)
S
r
1
r
2
=
M
∑
k
1
=
0
n
1
-
1
∑
k
2
=
0
n
2
-
1
exp
⁡
2
π
i
k
1
r
1
n
1
exp
⁡
2
π
i
k
2
r
2
n
2
X
k
1
k
2for fixed r1=0,1,…,n1-1 and r2=0,1,…,n2-1 discretizing κ1 and κ2. Squaring (10) yields Toeplitz matrices with Kronecker product structure (A.7) for the detected signal(11)Sr1r22=M2·xTr1⊗Tr2xwith multiple row and column indices (p1p2) and (q1q2), respectively, referring to the ordinary rows p1,p2 and columns q1,q2 of the Toeplitz matrices (9):(12)Tr1⊗Tr2p1p2,q1q2=Tr1p1q1Tr2p2q2,xq1q2≔Xq1q2.The scattering amplitudes in Fourier domain can be assigned with arbitrary phases. As this leaves the given intensity distributionS2 unchanged an related infinite set of broadened scattering data is formed by convolution. But if the signal S2 is assumed to be oversampled the original scattering data in this set appear to be sparse and can be recovered by an l1 minimization run on, for example, (8) or (11)—up to translations, reflections, or global phases.


## 3. Kalman Filter-Drivenl1 Minimization

### 3.1. Kalman Filter Equations

The equations for the Kalman filter usually apply to vectors over the field of real numbersR. We will examine next that the equations can also be extended to the field C. Anticipating this in the classical state space approach [36] the vectorial quantity xk is supposed to evolve according to the linear evolution (13a)xk+1=Akxk+uk+Gw∈Cnstochastic(13b)tochicxkxkxkxkAk∈Cn×n,uk∈Cndeterministic(13c)tochtixkxkxkxkxkkkxkG∈Cn×r,w∈Crstochasticwith a fixed matrix G and a determined sequence of evolution matrices Ak. In the model above the quantity xk is assumed to be only traceable indirectly by linear observations yk which can be viewed as linear mappings with given sensing matrices Ck by (14a)yk=Ckxk+v∈Cmstochastic(14b)CkCkstochasticxkxkCk∈Cm×n,v∈Cmstochastic.The stochastic behaviour of all the considered quantities is modelled by zero-mean Gaussian distributions with positive definite matricesR and Q. These covariance matrices and related mean values can be calculated from normalized Gaussian integrals denoted by expectation or mean values E∘:(15)Ew=0∈Cr,Covw=Eww=Q∈Cr×r,Ev=0∈Cm,Covv=Evv=R∈Cm×m.The main idea behind Kalman filtering is to invert the observation model (14a) and (14b), where the estimation of the quantities xk from the measurements yk shows predictor-corrector structure [36]: the prediction step relates the estimates xk+1- and xk+ by extrapolation(16a)xk+1-=Akxk++uk,Pk+=Covxk+,(16b)Pk+1-=AkPk+AkH+GQGH, whereas the correction step updates the estimate xk+1- by relating it to the new measurement yk+1: as the underlying equations utilize conditional mean values the update can be formulated in terms of covariances (17a)Pk+1+-1=Pk+1--1+Ck+1HR-1Ck+1,(17b)xk+1+=Pk+1+Pk+1--1xk+1-+Pk+1+Ck+1HR-1yk+1.An alternative form of the covariance cycle in the correction step above reads, in terms of the explicit so-called Kalman gain matrixKk+1, (18a)Kk+1=Pk+1-Ck+1HCk+1Pk+1-Ck+1H+R-1,(18b)Pk+1+=1-Kk+1Ck+1Pk+1-,(18c)xk+1+=xk+1-+Kk+1yk+1-Ck+1xk+1-.


### 3.2.l1 Minimization

The main problem in CS is to minimizex1 subject to the constraint y=Cx for fixed y and C. Because this linear structure matches the evolution equation (13a) of Kalman filtering in Section 3.1 we want to follow [29] using a pseudomeasurement to iteratively minimize the l1 norm: here, the kth estimate of the state vector x can be associated with xk and the norm minimization is driven by the fixed constraint y augmented by the lowered norm γxk-11 from the previous step by a factor of 0<γ≤1 as an additional observation. The main issue of applying the Kalman filter to the l1 minimization procedure consists in suitably linearizing the nonanalytic l1 norm: Solving(19)z2=Re2zz0¯z0+Im2zz0¯z0 for z in the vicinity of z0∈C yields for φ≔arg⁡(zz0¯/z0) the expression(20)z=Rezz0¯z01+tan2⁡φ≥Rezz0¯z0.Thus thel1 norm of any vector z∈Cn can be linearized in the vicinity of z0∈Cn utilizing its phase information p∈Cn by(21)Rep∣z≤z1,p=pz0≔z0¯1z01,…,z0¯nz0n.Due to vectors over the fieldC it is necessary to distinguish row vectors p from the usual column vectors p≔p related by Hermitian conjugation. The notation is based on complex scalar products and matrix multiplication (cf. Appendix A).


### 3.3. Pseudomeasurements

Involving the observationy=Cx∈Cm as a constraint like, for example, a Fourier transform, the Kalman filter equations (16a), (16b), (17a), and (17b) read (for A=G=1, uk=0) for the estimates xk≔xk+ to the vector x involving a full prediction-correction step (22a)Pk+1-1=Pk+Q-1+Ck+1HR-1Ck+1(22b)xk+1=Pk+1Pk-1xk+Pk+1Ck+1HR-1yk+1with the initial values P0∈Cn×n, x0∈Cn and the constant parameters C∈Cm×n, R∈C(m+1)×n and Q∈Cn×n to reconstruct x=limk→∞xk from a given y as a limiting value. The minimization of the l1 norm is incorporated into the (pseudo)measurements by(23)yk+1=yγkxk1∈Cm+1,Ck+1=Cpxk∈Cm+1×n, where the factor 0<γk≤1 adaptively lowers the current l1 norm solving for a corresponding estimate xk+1 in the next iteration step. For a solution to the linear constraint C∈Cm×n as initial data can serve(24)x0=CHC-1CHyform>nCHCCH-1yform≤n.


### 3.4. Convergence Considerations

By inserting (22a) into (22b) the enhancement of the estimate xk+1 compared to xk reads with respect to the measurement yk+1∈Cm+1 of a lowered l1 norm(25)Ck+1xk+1-yk+1=R-Ck+1Pk+1Ck+1HR-1·Ck+1xk-yk+1.Using the covariance cycle (22a) again the enhancement matrix composed of positive definite matrices R, P, and Q on the RHS can be recast into the more suitable3 form(26)R-Ck+1Pk+1Ck+1HR-1=R·R+Ck+1Pk+QCk+1H-1≤1m+1yielding the inequality(27)ej∣Ck+1xk+1-yk+1≤ej∣Ck+1xk-yk+1,j=1,…,m+1,by components. The key ingredients are the restricted positive eigenvalues of (26) bounded from above by the spectral norm 1m+12=1 (cf. (B.1)). What remains to prove for the inequality in (26) is the combination Ck(P+Q)CkH to be positive semidefinite.In the overdetermined casem≥n this holds true for CkHCk to be positive definite: as the combination Ck(P+Q)CkH∈C(m+1)×(m+1) has only rank n the amount of m+1-n of its eigenvalues is zero. The positive sign of the remaining eigenvalues can be obtained from a Cholesky factorization of P+Q=FFH with a lower triangular matrix F∈Cn×n yielding CkFFHCkH=(CkF)(CkF)H. Due to a singular value decomposition [37] of CkF the remaining eigenvalues belong to (CkF)H(CkF)=FH(CkHCk)F∈Cn×n which, however, is along with the prerequisite CkHCk>0 positive definite because of Theorem 1.Theorem 1 (see [37]).
IfB∈Cn×n is positive definite and X∈Cn×k has rank k, then the matrix Y=XHBX∈Ck×k is also positive definite.By virtue of the same theorem the combinationCk(P+Q)CkH∈C(m+1)×(m+1) is always positive definite in the underdetermined case m<n, which is independent of the sensing matrix Ck.So in each iteration step with adaptive0<γk≤1 the linearized l1 norm is lowered (or remains at least constant) subject to the approximated constraint y=Ckxk. In the real case (20) even holds with “=” yielding for j=m+1 the relation xk+11-γkxk1≤1-γkxk1 implying xk+11≤xk1 for all γk. To show a general convergence let us start with the exact solution x∞=x and C∞ related to it. Because of A=G=1 with zero shifts uk=0 in the Kalman filter model (13a), (13b), and (13c) the constant (pseudo)measurement y∞ already implies with γk≡1 a linear convergence of the series C∞x0,C∞x1,C∞x2,… by(28)ej∣C∞xk+1-y∞≤ej∣C∞xk-y∞≤⋯≤ej∣C∞x0-y∞for all j=1,…,m+1 starting with an x0 in the vicinity of x∞. Thus the series x0,x1,x2,… is said to be weakly convergent and reconstructs the original signal as long as the requirement for sparsity in the CS approach is met. In the framework of weak convergence the limiting value of Covxk reads for large k≫1(29)Pk≥P0-P0C∞HC∞P0C∞H-1C∞P0+O1k,m<n1kC∞HR-1C∞-1+O1k2,m≥n,where “=” represents the special case Q≡0. The limit Q→∞ is only possible for m≥n yielding a constant covariance Pk=(C∞HR-1C∞)-1 in each iteration.For an example of linear compressive sensing in the framework of Kalman filtering see [30, 31] using random samples of Fourier coefficients. Applying the method to a quadratic nonlinear observation model is shown for coherent X-ray scattering in Section 5.


## 3.1. Kalman Filter Equations

The equations for the Kalman filter usually apply to vectors over the field of real numbersR. We will examine next that the equations can also be extended to the field C. Anticipating this in the classical state space approach [36] the vectorial quantity xk is supposed to evolve according to the linear evolution (13a)xk+1=Akxk+uk+Gw∈Cnstochastic(13b)tochicxkxkxkxkAk∈Cn×n,uk∈Cndeterministic(13c)tochtixkxkxkxkxkkkxkG∈Cn×r,w∈Crstochasticwith a fixed matrix G and a determined sequence of evolution matrices Ak. In the model above the quantity xk is assumed to be only traceable indirectly by linear observations yk which can be viewed as linear mappings with given sensing matrices Ck by (14a)yk=Ckxk+v∈Cmstochastic(14b)CkCkstochasticxkxkCk∈Cm×n,v∈Cmstochastic.The stochastic behaviour of all the considered quantities is modelled by zero-mean Gaussian distributions with positive definite matricesR and Q. These covariance matrices and related mean values can be calculated from normalized Gaussian integrals denoted by expectation or mean values E∘:(15)Ew=0∈Cr,Covw=Eww=Q∈Cr×r,Ev=0∈Cm,Covv=Evv=R∈Cm×m.The main idea behind Kalman filtering is to invert the observation model (14a) and (14b), where the estimation of the quantities xk from the measurements yk shows predictor-corrector structure [36]: the prediction step relates the estimates xk+1- and xk+ by extrapolation(16a)xk+1-=Akxk++uk,Pk+=Covxk+,(16b)Pk+1-=AkPk+AkH+GQGH, whereas the correction step updates the estimate xk+1- by relating it to the new measurement yk+1: as the underlying equations utilize conditional mean values the update can be formulated in terms of covariances (17a)Pk+1+-1=Pk+1--1+Ck+1HR-1Ck+1,(17b)xk+1+=Pk+1+Pk+1--1xk+1-+Pk+1+Ck+1HR-1yk+1.An alternative form of the covariance cycle in the correction step above reads, in terms of the explicit so-called Kalman gain matrixKk+1, (18a)Kk+1=Pk+1-Ck+1HCk+1Pk+1-Ck+1H+R-1,(18b)Pk+1+=1-Kk+1Ck+1Pk+1-,(18c)xk+1+=xk+1-+Kk+1yk+1-Ck+1xk+1-.


## 3.2.l1 Minimization

The main problem in CS is to minimizex1 subject to the constraint y=Cx for fixed y and C. Because this linear structure matches the evolution equation (13a) of Kalman filtering in Section 3.1 we want to follow [29] using a pseudomeasurement to iteratively minimize the l1 norm: here, the kth estimate of the state vector x can be associated with xk and the norm minimization is driven by the fixed constraint y augmented by the lowered norm γxk-11 from the previous step by a factor of 0<γ≤1 as an additional observation. The main issue of applying the Kalman filter to the l1 minimization procedure consists in suitably linearizing the nonanalytic l1 norm: Solving(19)z2=Re2zz0¯z0+Im2zz0¯z0 for z in the vicinity of z0∈C yields for φ≔arg⁡(zz0¯/z0) the expression(20)z=Rezz0¯z01+tan2⁡φ≥Rezz0¯z0.Thus thel1 norm of any vector z∈Cn can be linearized in the vicinity of z0∈Cn utilizing its phase information p∈Cn by(21)Rep∣z≤z1,p=pz0≔z0¯1z01,…,z0¯nz0n.Due to vectors over the fieldC it is necessary to distinguish row vectors p from the usual column vectors p≔p related by Hermitian conjugation. The notation is based on complex scalar products and matrix multiplication (cf. Appendix A).


## 3.3. Pseudomeasurements

Involving the observationy=Cx∈Cm as a constraint like, for example, a Fourier transform, the Kalman filter equations (16a), (16b), (17a), and (17b) read (for A=G=1, uk=0) for the estimates xk≔xk+ to the vector x involving a full prediction-correction step (22a)Pk+1-1=Pk+Q-1+Ck+1HR-1Ck+1(22b)xk+1=Pk+1Pk-1xk+Pk+1Ck+1HR-1yk+1with the initial values P0∈Cn×n, x0∈Cn and the constant parameters C∈Cm×n, R∈C(m+1)×n and Q∈Cn×n to reconstruct x=limk→∞xk from a given y as a limiting value. The minimization of the l1 norm is incorporated into the (pseudo)measurements by(23)yk+1=yγkxk1∈Cm+1,Ck+1=Cpxk∈Cm+1×n, where the factor 0<γk≤1 adaptively lowers the current l1 norm solving for a corresponding estimate xk+1 in the next iteration step. For a solution to the linear constraint C∈Cm×n as initial data can serve(24)x0=CHC-1CHyform>nCHCCH-1yform≤n.


## 3.4. Convergence Considerations

By inserting (22a) into (22b) the enhancement of the estimate xk+1 compared to xk reads with respect to the measurement yk+1∈Cm+1 of a lowered l1 norm(25)Ck+1xk+1-yk+1=R-Ck+1Pk+1Ck+1HR-1·Ck+1xk-yk+1.Using the covariance cycle (22a) again the enhancement matrix composed of positive definite matrices R, P, and Q on the RHS can be recast into the more suitable3 form(26)R-Ck+1Pk+1Ck+1HR-1=R·R+Ck+1Pk+QCk+1H-1≤1m+1yielding the inequality(27)ej∣Ck+1xk+1-yk+1≤ej∣Ck+1xk-yk+1,j=1,…,m+1,by components. The key ingredients are the restricted positive eigenvalues of (26) bounded from above by the spectral norm 1m+12=1 (cf. (B.1)). What remains to prove for the inequality in (26) is the combination Ck(P+Q)CkH to be positive semidefinite.In the overdetermined casem≥n this holds true for CkHCk to be positive definite: as the combination Ck(P+Q)CkH∈C(m+1)×(m+1) has only rank n the amount of m+1-n of its eigenvalues is zero. The positive sign of the remaining eigenvalues can be obtained from a Cholesky factorization of P+Q=FFH with a lower triangular matrix F∈Cn×n yielding CkFFHCkH=(CkF)(CkF)H. Due to a singular value decomposition [37] of CkF the remaining eigenvalues belong to (CkF)H(CkF)=FH(CkHCk)F∈Cn×n which, however, is along with the prerequisite CkHCk>0 positive definite because of Theorem 1.Theorem 1 (see [37]).
IfB∈Cn×n is positive definite and X∈Cn×k has rank k, then the matrix Y=XHBX∈Ck×k is also positive definite.By virtue of the same theorem the combinationCk(P+Q)CkH∈C(m+1)×(m+1) is always positive definite in the underdetermined case m<n, which is independent of the sensing matrix Ck.So in each iteration step with adaptive0<γk≤1 the linearized l1 norm is lowered (or remains at least constant) subject to the approximated constraint y=Ckxk. In the real case (20) even holds with “=” yielding for j=m+1 the relation xk+11-γkxk1≤1-γkxk1 implying xk+11≤xk1 for all γk. To show a general convergence let us start with the exact solution x∞=x and C∞ related to it. Because of A=G=1 with zero shifts uk=0 in the Kalman filter model (13a), (13b), and (13c) the constant (pseudo)measurement y∞ already implies with γk≡1 a linear convergence of the series C∞x0,C∞x1,C∞x2,… by(28)ej∣C∞xk+1-y∞≤ej∣C∞xk-y∞≤⋯≤ej∣C∞x0-y∞for all j=1,…,m+1 starting with an x0 in the vicinity of x∞. Thus the series x0,x1,x2,… is said to be weakly convergent and reconstructs the original signal as long as the requirement for sparsity in the CS approach is met. In the framework of weak convergence the limiting value of Covxk reads for large k≫1(29)Pk≥P0-P0C∞HC∞P0C∞H-1C∞P0+O1k,m<n1kC∞HR-1C∞-1+O1k2,m≥n,where “=” represents the special case Q≡0. The limit Q→∞ is only possible for m≥n yielding a constant covariance Pk=(C∞HR-1C∞)-1 in each iteration.For an example of linear compressive sensing in the framework of Kalman filtering see [30, 31] using random samples of Fourier coefficients. Applying the method to a quadratic nonlinear observation model is shown for coherent X-ray scattering in Section 5.


## 4. Linearized Observation Model

Like thel1 norm in Section 3.2 we linearize the squared observations (8) to apply the Kalman filtering scheme. Differentiating with respect to xk and xk¯, respectively, yields(30)∂∂xkSr2=xTrk,∂∂xk¯Sr2=Trxkand Taylor expansion around z0∈Cn reads up to 1st order(31)Sr2=2Rez0Trx-z0Trz0+⋯,r=0,1,…,n-1.With the biasessr=z0Trz0,r=0,1,…,n-1, from the linearization the observation model is(32)S02⋮Sn-12x1=2z0T0⋮2z0Tn-1pz0·x0⋮xn-2xn-1-s0⋮sn-10(33)yk=Ck·x-sk,where the real parts over all Cx combinations have to be taken. The model in detail reads yk=Re(Ckx)-sk∈Rn+1 and because of the Cx combinations we choose the alternative representation (18a), (18b), and (18c) of the Kalman filter equations to invert (32). Then the iterated estimates xk to the state vector x read (34a)Kk+1=Pk+QCk+1HCk+1Pk+QCk+1H+R-1,(34b)Pk+1=1-Kk+1Ck+1Pk+Q,(34c)xk+1=xk+Kk+1·Reyk+1-Ck+1xk+sk+1with the complex sensing matrix Ck+1∈C(n+1)×n, biases sk+1∈Cn+1, and the given real (pseudo)measurements yk+1∈Rn+1 according to(35)Ck+1=2xkT0⋮2xkTn-1pxk,sk+1=xkT0xk⋮xkTn-1xk0,yk+1=S02⋮Sn-12γkxk1.A properly chosen factor0<γk≤1 adaptively lowers the l1 norm in each iteration step k to reconstruct the signal x=limk→∞xk as a limiting value from its given squared measurements y=(S02,…,Sn-12). Note that T0xk,T1xk,…,Tn-1xk⊂Cn form an orthogonal basis. Thus(36)Ck+1HCk+1=4∑r=0n-1TrxkxkTr+ppxkis Hermitian and positive definite due to (B.1) for all k with xk≠0 which is sufficient to prove weak convergence of the l1 minimization under the constraint of squared observations (cf. (26)–(28)).


### 4.1. Example 1: Stacking Sequence

To demonstrate the relative phase recovery from the observed intensity (8) not only restricting to the Zinc-Blende and Wurtzite phases (6) assume a linear chain with 117 lattice sites and 6 equidistant sparse amplitudes xj=1 with phases from the set -π/3,-π/6,0,π/6,π/3,π/2 forming xsparse∈Cn.As in the original objective the number of bilayers is roughly known we use as initial valuesx0∈Cn a broadened modulus and phase distribution related to the setting xsparse (cf. Figure 3): with the support S≔supp⁡xsparse a leakage(37)f0±k=∑α∈Ssin⁡k±s-απk±s-απ∈R,k=0,…,n-1for the moduli and phases of xsparse is with shifts 0≤s<1, a suitable way to also consider the limit s→0 of exactly known positions. Note that for the phases we only used one fixed value for each occupied lattice site broadened by the leakage. So the main effort is the recovery of the relative phases from the measurements.Figure 3
A number of6 equidistant scatterers on a periodic chain of length 117. To allow for a signal from the substrate the wires are growing on we added some parasitic scatterers to the chain. As initial data we use the moduli and phases from xsparse symmetrically broadened by a leakage 1/2(f0++f0-) (cf. (37)) with shift s=0.3 and amplitudes 0.6 and 1.1, respectively.For an exponential decaying lowering factorγk=1-0.1·exp⁡(-0.0019k) within 1200 iterations the reconstruction results for noiseless synthetic measurements are shown in Figures 4 and 5. Experientially, good algorithm’s covariances were found to be(38)P0=0.3·1117,Q=10-8·1117,R=diag⁡10-4,…,10-4,10-6∈R118×118.Figure 4
Reconstructed moduli of all scatterers. In the phase domain the values corresponding to vanishing moduli are suppressed. Due to the symmetry of the sensing problem the phases are only reconstructed up to a global phase.Figure 5
Adding a global phase (slider on the left) verifies that all relative phases are reconstructed correctly. The open rectangles mark the position of occupied lattice sites on the chain. In the original objective this would refer to bilayers and substrate in the nanowire.To drive thel1 minimization the corresponding entry (here 10-6) in the R matrix has to be much lower in magnitude than the ones (here 10-4) related to the observations of the constraint. When reconstructed amplitudes are said to vanish (e.g., in Figure 4) this means only up to numerical precision where the zero value is dominated by the entries of Q. Empirically, the exponential decay in the lowering factor γk is chosen such that its values approach about 0.99,…,0.999 after the maximum number of iterations. A resulting typical smooth convergence of the l1 norm due to the exponential decay is shown in Figure 6: the constant tail can be used to determine a maximal number of iterations as a stopping criterion. Note that the reconstructed results beyond this guessed number are insensitive to further iterations meaning for (34a), (34b), and (34c) to perform fixed point iterations within accuracies related to the covariances Q and Pk.Figure 6
Note that thel1 convergence (27) in the kth iteration step relates the linearized norm p∣xk≤xk1 to xk-11. Because of this inequality (cf. (21)) there is no guarantee of a monotonically lowered l1 norm, which can be seen in the beginning (iteration step 2).In the nanowire the equidistant bilayers are at fixed positions in growth directiona3. Thus the reconstruction of their complex scattering amplitudes does not suffer from “off-grid” problems in general (see, e.g., [38–40]) if integer multiples of the lattice constants are used for the DFT.


### 4.2. Example 2: Pattern Reconstruction

The 2D reconstruction from (11) can be calculated with the same algorithm already used for Example 1 above, as the 2D data set of size 12×25 can be mapped to a 1D vector by means of (12) (cf. Figure 7). Like the example from Figure 3, we use as initial values with the same reasoning and broadening parameters the given distribution. Within 2000 iterations and an exponentially decaying lowering factor γk=1-0.17·exp⁡(-0.0028k) the reconstruction for noiseless measurements is shown in Figure 8 with empirically good algorithm’s covariances:(39)P0=0.3·1300,Q=10-7·1300,R=diag⁡10-4,…,10-4,10-6∈R301×301.Note that there is a combined translational and reflectional symmetry mapping the vectorized real scattering amplitude distribution onto itself (cf. Figure 7). As this is also a symmetry of the squared observations (8) in the linear configuration, there are two competing solutions to the sensing problem. Interestingly, breaking this symmetry improved the l1 convergence in Figure 8: without the parasitic scatterers the symmetric positions of consecutive sites4 are correctly (with respect to the support) found after 5000 iterations (cf. Figure 9) and the remaining discrepancies in the amplitudes do not even vanish completely after inefficient 50000 iterations. On the contrary the random sample of 25 scatterers from Figure 10 is already reconstructed with less than 1000 iterations suggesting a strong dependence of the algorithm’s convergence from the vector under consideration.Figure 7
The actual scatterers forming the5×5 structure are made from real amplitudes 0.12, 0.28, and 0.56 whereas the two parasitic scatterers are assigned with 0.23 and 0.42 (a). The scattering strength in the plane is represented by the area of the filled circles (b).
(a) (b)Figure 8
The real amplitudes are nicely reconstructed (a) after2000 iterations which can also be seen in the 2D setting (b) by comparing with Figure 7.
(a) (b)Figure 9
Without parasitic scatterers the positions are correctly found (b) after5000 iterations and γk=1-0.17⋅exp⁡(-0.0012k). The remaining discrepancies in the amplitudes seem to be complementarily symmetric with respect to the central peak (a).
(a) (b)Figure 10
A random distribution of25 scatterers with real random amplitudes normalized in the l2 norm. The reconstruction is without artefacts and discrepancies occurred after 1000 iterations with γk=1-0.17⋅exp⁡(-0.0058k).


## 4.1. Example 1: Stacking Sequence

To demonstrate the relative phase recovery from the observed intensity (8) not only restricting to the Zinc-Blende and Wurtzite phases (6) assume a linear chain with 117 lattice sites and 6 equidistant sparse amplitudes xj=1 with phases from the set -π/3,-π/6,0,π/6,π/3,π/2 forming xsparse∈Cn.As in the original objective the number of bilayers is roughly known we use as initial valuesx0∈Cn a broadened modulus and phase distribution related to the setting xsparse (cf. Figure 3): with the support S≔supp⁡xsparse a leakage(37)f0±k=∑α∈Ssin⁡k±s-απk±s-απ∈R,k=0,…,n-1for the moduli and phases of xsparse is with shifts 0≤s<1, a suitable way to also consider the limit s→0 of exactly known positions. Note that for the phases we only used one fixed value for each occupied lattice site broadened by the leakage. So the main effort is the recovery of the relative phases from the measurements.Figure 3
A number of6 equidistant scatterers on a periodic chain of length 117. To allow for a signal from the substrate the wires are growing on we added some parasitic scatterers to the chain. As initial data we use the moduli and phases from xsparse symmetrically broadened by a leakage 1/2(f0++f0-) (cf. (37)) with shift s=0.3 and amplitudes 0.6 and 1.1, respectively.For an exponential decaying lowering factorγk=1-0.1·exp⁡(-0.0019k) within 1200 iterations the reconstruction results for noiseless synthetic measurements are shown in Figures 4 and 5. Experientially, good algorithm’s covariances were found to be(38)P0=0.3·1117,Q=10-8·1117,R=diag⁡10-4,…,10-4,10-6∈R118×118.Figure 4
Reconstructed moduli of all scatterers. In the phase domain the values corresponding to vanishing moduli are suppressed. Due to the symmetry of the sensing problem the phases are only reconstructed up to a global phase.Figure 5
Adding a global phase (slider on the left) verifies that all relative phases are reconstructed correctly. The open rectangles mark the position of occupied lattice sites on the chain. In the original objective this would refer to bilayers and substrate in the nanowire.To drive thel1 minimization the corresponding entry (here 10-6) in the R matrix has to be much lower in magnitude than the ones (here 10-4) related to the observations of the constraint. When reconstructed amplitudes are said to vanish (e.g., in Figure 4) this means only up to numerical precision where the zero value is dominated by the entries of Q. Empirically, the exponential decay in the lowering factor γk is chosen such that its values approach about 0.99,…,0.999 after the maximum number of iterations. A resulting typical smooth convergence of the l1 norm due to the exponential decay is shown in Figure 6: the constant tail can be used to determine a maximal number of iterations as a stopping criterion. Note that the reconstructed results beyond this guessed number are insensitive to further iterations meaning for (34a), (34b), and (34c) to perform fixed point iterations within accuracies related to the covariances Q and Pk.Figure 6
Note that thel1 convergence (27) in the kth iteration step relates the linearized norm p∣xk≤xk1 to xk-11. Because of this inequality (cf. (21)) there is no guarantee of a monotonically lowered l1 norm, which can be seen in the beginning (iteration step 2).In the nanowire the equidistant bilayers are at fixed positions in growth directiona3. Thus the reconstruction of their complex scattering amplitudes does not suffer from “off-grid” problems in general (see, e.g., [38–40]) if integer multiples of the lattice constants are used for the DFT.


## 4.2. Example 2: Pattern Reconstruction

The 2D reconstruction from (11) can be calculated with the same algorithm already used for Example 1 above, as the 2D data set of size 12×25 can be mapped to a 1D vector by means of (12) (cf. Figure 7). Like the example from Figure 3, we use as initial values with the same reasoning and broadening parameters the given distribution. Within 2000 iterations and an exponentially decaying lowering factor γk=1-0.17·exp⁡(-0.0028k) the reconstruction for noiseless measurements is shown in Figure 8 with empirically good algorithm’s covariances:(39)P0=0.3·1300,Q=10-7·1300,R=diag⁡10-4,…,10-4,10-6∈R301×301.Note that there is a combined translational and reflectional symmetry mapping the vectorized real scattering amplitude distribution onto itself (cf. Figure 7). As this is also a symmetry of the squared observations (8) in the linear configuration, there are two competing solutions to the sensing problem. Interestingly, breaking this symmetry improved the l1 convergence in Figure 8: without the parasitic scatterers the symmetric positions of consecutive sites4 are correctly (with respect to the support) found after 5000 iterations (cf. Figure 9) and the remaining discrepancies in the amplitudes do not even vanish completely after inefficient 50000 iterations. On the contrary the random sample of 25 scatterers from Figure 10 is already reconstructed with less than 1000 iterations suggesting a strong dependence of the algorithm’s convergence from the vector under consideration.Figure 7
The actual scatterers forming the5×5 structure are made from real amplitudes 0.12, 0.28, and 0.56 whereas the two parasitic scatterers are assigned with 0.23 and 0.42 (a). The scattering strength in the plane is represented by the area of the filled circles (b).
(a) (b)Figure 8
The real amplitudes are nicely reconstructed (a) after2000 iterations which can also be seen in the 2D setting (b) by comparing with Figure 7.
(a) (b)Figure 9
Without parasitic scatterers the positions are correctly found (b) after5000 iterations and γk=1-0.17⋅exp⁡(-0.0012k). The remaining discrepancies in the amplitudes seem to be complementarily symmetric with respect to the central peak (a).
(a) (b)Figure 10
A random distribution of25 scatterers with real random amplitudes normalized in the l2 norm. The reconstruction is without artefacts and discrepancies occurred after 1000 iterations with γk=1-0.17⋅exp⁡(-0.0058k).


## 5. Conclusion and Outlook

We aimed to implement nonlinear compressive sensing for complex vectorsx by inverting the underlying observation model with Kalman filtering. For this reason we proved a weak convergence of the filter equations for complex sensing matrices C and Hermitian covariances R,P,Q. For the example of quadratic nonlinearities we applied our formulas to retrieve relative phase information in the objective of simulated noiseless coherent X-ray diffraction. Due to the nonlinearity the noise in the intensity is not Gaussian and needs further considerations to also account for, for example, lattice distortions or the finite coherence of the primary beam. As an outlook we presented a 2D pattern reconstruction which, hopefully, can help to investigate if the cross sections of the wires were grown regularly.Because of the Jacobians building up the sensing matrixC the convergence of the underlying l1 norm minimization does depend on the reconstructed vector x and may go beyond the resolution issues [40] for even constant sensing matrices. For this reason a thorough analysis on the relation between maximal sparsity (for the examples here about 10% of the available lattice sites) and the Toeplitz matrix (9) representing the sensing process for a successful CS is needed. As mentioned in the introduction algorithms in general minimizing the l1 norm could be of interest retrieving information out from intensity spectra: to compare our results we applied the nonlinear version [41] of the primal dual algorithm [42] to the 1D sensing problem above yielding similar results with respect to iteration numbers, reconstructed amplitudes, and phases. On the one hand primal dual reconstructed the zeros outside the support (once a solution was isolated from the algorithm) numerically exact compared to accuracies of orders 10-5 reached by the Kalman filter. On the other hand our approach focused in the first quarter of the iterations on guessing the support by strongly increasing and decreasing the corresponding amplitudes, whereas primal dual homogeneously acted on all the amplitudes during all iterations. For a detailed comparison and maybe a combination, the whole parameter range of the approaches have to be investigated. To apply the Kalman filter-based algorithm above to recorded data [17] we need to deal with several 104x-vector’s entries describing a sparsely (by a factor of roughly 1 : 101) occupied linear grid covering all the approximately 103 bilayers in a nanowire of 500 nm in height. For this reason matrix inversions, cf. (34a), (34b), and (34c), should be reformulated by, for example, sequential processing techniques to make the algorithm more scalable. Furthermore it would also be interesting if other adaptive lowering factors γk compared to the exponential ones yield faster convergences.As the reconstruction from quadratic constraints seems to depend only little on the explicit algorithm used for thel1 minimization, it would be interesting to figure out to which extent the linear matrix completion approach [22–24] with respect to the nuclear norm can still be applied in the face of seemingly too few numbers of independent observations.



---

*Source: 9853714-2016-07-05.xml*