Semiparametric estimation for incoherent optical imaging

Semiparametric estimation f or incoherent optical imaging Mankei Tsang ∗ Department of Electrical and Computer Engineering, National University of Singapore , 4 Engineering Drive 3, Singapor e 117583 and Department of Physics, National University of Singapor e, 2 Science Drive 3, Singapor e 117551 (Dated: June 18, 2019) The theory of semiparametric estimation offers an elegant way of computing the Cram ´ er-Rao bound for a parameter of interest in the midst of infinitely many nuisance parameters. Here I apply the theory to the problem of moment estimation for incoherent imaging under the effects of diffraction and photon shot noise. Using a Hilbert-space formalism designed for Poisson processes, I deri ve e xact semiparametric Cram ´ er-Rao bounds and efficient estimators for both direct imaging and a quantum-inspired measurement method called spatial-mode demultiplexing (SP ADE). The results establish the superiority of SP ADE ev en when little prior information about the object is av ailable. I. INTR ODUCTION T w o fundamental problems confront incoherent optical imaging: the diffraction limit [ 1 , 2 ] and the photon shot noise [ 3 , 4 ]. T o quantify their ef fects on the resolution rigorously , the Cram ´ er-Rao bound (CRB) on the error of parameter es- timation [ 5 ] has been widely used, especially in astronomy and fluorescence microscopy [ 6 – 17 ]. Most previous studies, howe ver , assume that the object has a simple specific shape, such as a point source or two, and only one or few param- eters of the object are unkno wn. Such parametric models may not be justifiable when there is little prior information about the object. W ithout a parametric model, the CRB seems intractable—infinitely many parameters are needed to specify the object distribution, leading to a Fisher information matrix with infinitely many entries, and then the infinite-dimensional matrix has to be in verted to giv e the CRB. While there also exist man y studies on superresolution that can deal with more general objects [ 17 – 20 ], they either ignore noise or use noise models that are too simplistic to capture the signal-dependent nature of photon shot noise. T o compute the CRB and to ev aluate the efficiency of es- timators for general objects, here I propose a theory of semi- parametric estimation for incoherent optical imaging. Semi- parametric estimation refers to the estimation of a parameter of interest in the presence of infinitely many other unknown “nuisance” parameters [ 21 , 22 ]. The method has found many applications in econometrics, biostatistics, and astrostatistics [ 21 ]. A typical example is the estimation of the mean of a ran- dom variable when its probability density is assumed to ha ve finite variance but otherwise arbitrary . Thanks to a beautiful Hilbert-space formalism [ 21 , 22 ], the semiparametric theory is able to compute the CRB for such problems despite the in- finite dimensionality and also e v aluate the existence and effi- ciency of semiparametric estimators. Such problems are ex- actly the type that bedevil the study of imaging thus far , and here I show ho w the semiparametric theory can be used to yield similarly elegant results for optical imaging. The optics problem of interest here is the far-field imag- ing of an object emitting spatially incoherent light [ 2 , 4 ], with ∗ mankei@nus.edu.sg ; https://www .ece.nus.edu.sg/stfpage/tmk/ the most important applications being optical astronomy [ 6 – 9 ] and fluoresence microscopy [ 10 – 14 ]. W ith a finite numer- ical aperture, the imaging system introduces a spatial band- width limit to the wav es, otherwise known as the dif fraction limit [ 1 , 2 ]. The standard measurement, called direct imag- ing, records the intensity of the light on the image plane. Re- cently , quantum information theory inspired the in vention of an alternativ e measurement called spatial-mode demultiplex- ing (SP ADE) [ 23 ], which has been shown theoretically [ 23 – 51 ] and experimentally [ 52 – 58 ] to be superior to direct imag- ing in resolving two sub-Rayleigh sources and estimating the size and moments of a subdiffraction object. Most of the aforementioned studies, ho wever , assume parametric mod- els for the object. Exceptions include Refs. [ 24 – 27 , 50 , 55 ], which consider the estimation of object moments, but the re- sults there are not conclusiv e—only the CRB for direct imag- ing was computed exactly [ 25 ], while the CRB for SP ADE was ev aluated only approximately [ 24 , 25 , 27 ]. Another prob- lem is the existence and efficiency of unbiased moment esti- mators; again only approximate results ha ve been obtained so far [ 24 , 25 ]. Building on the established semiparametric the- ory [ 21 , 22 ], here I compute the exact semiparametric CRBs and also propose unbiased and efficient moment estimators for both direct imaging and SP ADE. These results enable a fair and rigorous comparison of the two measurement meth- ods, which proves the fundamental superiority of SP ADE for moment estimation. This paper is organized as follows. Section II introduces the Fisher information and the CRB for Poisson processes. Section III presents the semiparametric CRB in terms of a Hilbert-space formalism designed for such processes. Sec- tion IV introduces the models of direct imaging and SP ADE. Section V computes the CRB for moment estimation with di- rect imaging and proposes an efficient estimator . Section VI shows ho w the CRB should be modified for a normalized ob- ject distribution. Section VII computes the CRB for SP ADE and also proposes an ef ficient estimator . Section VIII uses the CRBs to compare the performances of direct imaging and SP ADE, demonstrating the superiority of SP ADE for subd- iffraction objects. Section IX concludes the paper and points out open issues, while the Appendices detail the technical is- sues that arise in the main text. 2 II. CRAM ´ ER-RA O BOUND FOR POISSON PROCESSES For optical astronomy [ 4 , 6 , 8 , 9 ], fluorescence microscopy [ 10 – 14 ], and ev en electron microscopy [ 15 , 16 ], Poisson noise can be safely assumed. Suppose that each detector in a pho- todetector array is labeled by x ∈ X , where X denotes the detector space. Assume that the observed process, such as the image recorded by a camera, is a Poisson random mea- sure n on X and its σ -algebra Σ , with a mean given by the intensity measure ¯ n on the same ( X , Σ) [ 59 ]. n ( A ) for any A ∈ Σ is then a Poisson v ariable with mean ¯ n ( A ) . For ex- ample, if X ⊆ R 2 is a two-dimensional surface, then Σ is the set of all subareas that can be defined on the surface, and n ( A ) = R x ∈A dn ( x ) is the detected photon number over the area A . For an y v ectoral function h : X → R q on the detector space, ˇ h ( n ) = Z h ( x ) dn ( x ) , (2.1) a linear functional of n , is a random variable with statistics E ( ˇ h ) = Z h ( x ) d ¯ n ( x ) = ν ( h ) , (2.2) V ( ˇ h ) = E ( ˇ h ˇ h > ) − E ( ˇ h ) E ( ˇ h > ) = ν ( hh > ) , (2.3) where E denotes the statistical e xpectation, V denotes the co- variance, ν denotes the average with respect to the intensity measure ¯ n , > denotes the matrix transpose, and all vectors in this paper are column vectors. Suppose that ¯ n depends on an unknown vectoral parame- ter θ ∈ Θ ⊂ R p with p entries and has a density f ( x | θ ) with respect to a dominating measure µ such that f ( x | θ ) = d ¯ n ( x | θ ) /dµ ( x ) . The log-likelihood deriv ativ es are given by [ 60 ] ˇ S j ( n | θ ) = Z ∂ ∂ θ j ln f ( x | θ ) dn ( x ) − ∂ ∂ θ j Z d ¯ n ( x | θ ) . (2.4) As ˇ S is a linear functional of n , its cov ariance, called the Fisher information matrix, can be simplified via Eq. ( 2.3 ) and is giv en by [ 60 ] J = V ( ˇ S ) = Z S ( x )[ S ( x )] > d ¯ n ( x ) = ν ( S S > ) , (2.5) where S is a vector of detector-space functions gi ven by S j ( x | θ ) = ∂ ∂ θ j ln f ( x | θ ) . (2.6) Here V , ˇ S , ¯ n , S , and ν are all ev aluated at the same θ , and I assume hereafter that all functions of θ are ev aluated implic- itly at the same θ . Each S j is hereafter called a score function, borrowing the same terminology for ˇ S in statistics [ 21 , 22 ]. An important distinction is that, whereas E ( ˇ S j ) = 0 , ν ( S j ) does hav e to be zero, since ¯ n does not hav e to be normalized. Let β ( θ ) be a scalar parameter of interest. If β ( θ ) = θ k for example, then all the other parameters in θ are called nuisance parameters. For any unbiased estimator ˇ β ( n ) , the CRB on its variance is [ 5 ] V ( ˇ β ) ≥ u > J − 1 u = CRB , u j = ∂ β ∂ θ j . (2.7) J − 1 seems intractable if θ is infinite-dimensional. The next section introduces a clev erer method. III. SEMIP ARAMETRIC CRAM ´ ER-RA O BOUND The key to the semiparametric theory is to treat random variables as elements in a Hilbert space [ 21 , 22 ]. Here I in- troduce another Hilbert space for detector-space functions on top of the statistical one for the purpose of computing the CRB for Poisson processes. Define an inner product between two scalar functions h 1 , h 2 : X → R as h h 1 , h 2 i = ν ( h 1 h 2 ) = Z h 1 ( x ) h 2 ( x ) d ¯ n ( x ) , (3.1) and the norm as || h || = p h h, h i = p ν ( h 2 ) . (3.2) W ith the inner product, a Hilbert space H can be defined as the set of all square-summable functions, viz., H =  h ( x ) : ν ( h 2 ) < ∞  . (3.3) Denote the set of score functions { S j } as S in a slight abuse of notation. If the Fisher information J j j = ν ( S 2 j ) < ∞ for all j , S ⊂ H . Define the tangent space T ⊆ H of a parametric model as the linear span of S , or T =  w > S : w ∈ R p  = span( S ) . (3.4) Define also an “influence” function as any ˜ β ∈ H that satisfies ν ( ˜ β S ) = u, (3.5) borrowing the name of a similar concept in statistics [ 21 , 22 ]. The Cauchy-Schw artz inequality ν ( ˜ β 2 )[ w > ν ( S S > ) w ] ≥ ( u > w ) 2 with w = [ ν ( S S > )] − 1 u then yields ν ( ˜ β 2 ) ≥ u > J − 1 u, (3.6) the right-hand side of which coincides with the CRB giv en by Eq. ( 2.7 ). Define the ef ficient influence as the influence function that saturates Eq. ( 3.6 ), viz., ˜ β eff = u > J − 1 S = ν ( ˜ β S > )  ν ( S S > )  − 1 S, (3.7) CRB = ν ( ˜ β 2 eff ) . (3.8) Equation ( 3.7 ) can be interpreted as the orthogonal projection of any influence function ˜ β ∈ H that satisfies Eq. ( 3.5 ) into T , viz., ˜ β eff = Π( ˜ β |T ) = arg min h ∈T || ˜ β − h || . (3.9) 3 FIG. 1. The efficient influence ˜ β eff is the orthogonal projection of any influence function ˜ β ∈ H that satisfies Eq. ( 3.5 ) into the tangent space T = span( S ) . The norm of ˜ β eff giv es the CRB. Figure 1 illustrates this concept. Consider now the semiparametric scenario. For the purpose of this paper , it suf fices to assume that the dimension of θ is infinite but countable ( p = ∞ ). The score functions are still defined in the same way , but now there are infinitely many of them. The tangent space should be modified to be the closed linear span T = span( S ) , (3.10) so that projection into it is well defined [ 61 ], and the semi- parametric CRB is still given by Eqs. ( 3.5 ), ( 3.8 ), and ( 3.9 ); see Appendix A for a proof. This Hilbert-space approach is tractable when finding a candidate influence function is straightforward and the tangent space is so lar ge that the can- didate is already in it or at least very close to it. If the di- mension of θ is uncountable, the tangent space and the CRB can still be defined via the concept of parametric submodels [ 21 , 22 ], although it is not needed here. If β can be expressed as a functional β ( f ) , a useful way of finding an influence function is to consider a functional deriv ativ e of β ( f ) with respect to h ( x ) defined as ˙ β ( f , h ) = lim  → 0 β ((1 + h ) f ) − β ( f )  (3.11) = Z ˜ β ( x ) h ( x ) f ( x ) dµ ( x ) = ν ( ˜ β h ) , (3.12) which leads to ∂ β ∂ θ j = lim  → 0 β ( f + ∂ f /∂ θ j ) − β ( f )  (3.13) = ˙ β ( f , S j ) = ν ( ˜ β S j ) = u j , (3.14) and the ˜ β ( x ) function obtained from the functional deri v ativ e is an influence function that satisfies Eq. ( 3.5 ). The simplest example is the linear functional β ( f ) = Z ˜ β ( x ) f ( x ) dµ ( x ) = ν ( ˜ β ) , (3.15) and ˜ β ( x ) is an influence function. If the tangent space is so large that T = H , then a square- summable influence function is already in H = T and there- fore efficient. There are often some constraints that make T smaller , howe ver , and the CRB is reduced as a result. For example, if the constraint can be e xpressed as γ ( f ) = 0 , (3.16) and its functional deriv ativ e is ˙ γ ( f , h ) = ν ( h ˜ γ ) , (3.17) then ∂ γ ( f ) ∂ θ j = ˙ γ ( f , S j ) = ν ( ˜ γ S j ) = h ˜ γ , S j i = 0 , (3.18) and it follo ws that ˜ γ should be placed in the set that spans T ⊥ , the orthocomplement of T in H . In terms of T ⊥ , the ef ficient influence can be ev aluated as ˜ β eff = ˜ β − Π( ˜ β |T ⊥ ) . (3.19) If T ⊥ = span( R ) , then Π( ˜ β |T ⊥ ) = ν ( ˜ β R > )  ν ( RR > )  − 1 R, (3.20) which is still tractable if R has a low dimension. IV . INCOHERENT OPTICAL IMA GING Consider a distribution of spatially incoherent sources de- scribed by the measure F on the object plane with coordi- nate y , a far -field paraxial imaging system with point-spread function ψ ( z − y ) for the field [ 2 ], further processing of the field on the image plane with coordinate z via passiv e lin- ear optics with Green’ s function κ ( x, z ) , and Poisson noise at the output detectors labeled by x ∈ X , as depicted by Fig. 2 . For simplicity , assume one-dimensional imaging such that y , z ∈ R ; generalization for two-dimensional imaging is possible [ 24 , 25 ] but not very interesting. The intensity can be described by the mixture model [ 4 , 23 – 25 , 30 ] f ( x ) = Z     Z κ ( x, z ) ψ ( z − y ) dz     2 dF ( y ) , (4.1) where the image-plane coordinate z is normalized with re- spect to the magnification factor , both y and z are normal- ized with respect to the width of the point-spread function such that they are dimensionless, and ψ is normalized as R | ψ ( x ) | 2 dx = 1 . This semiclassical Poisson model can be deriv ed from standard quantum optics [ 23 , 51 , 62 ]. For direct imaging with infinitesimal pixels, κ ( x, z ) = √ τ δ ( x − z ) , where τ is a positi ve con version factor , x ∈ X = R denotes the position of each pixel, dµ ( x ) = dx , and the image intensity obeys the con volution model f ( x ) = Z H ( x − y ) dF ( y ) , H ( x ) = τ | ψ ( x ) | 2 , (4.2) 4 im age plane spatially incoher ent sources passiv e linea r optics detectors estim ator FIG. 2. A far-field incoherent imaging system. See the main text for definitions. which will be studied in Sec. V . The most remarkable physics of the problem lies in the possibility of improving the measurement via optics with a different Green’ s function κ . Quantum information theory has sho wn that substantial improv ement is possible for subd- iffraction objects, and SP ADE has been found to be quantum- optimal in many special cases [ 23 – 51 ]. In one v ersion of SP ADE, κ ∗ ( q , z ) is the q th mode function in the point-spread- function-adapted (P AD) basis [ 25 , 36 ], such that the output intensity is giv en by f ( q ) = Z H ( q | y ) dF ( y ) , q ∈ X = N 0 , (4.3) H ( q | y ) =     Z κ ( q , z ) ψ ( z − y ) dz     2 , (4.4) where µ is simply the counting measure and κ and H obey special properties, as further discussed in Sec. VII . For a fair comparison, the quantum efficiencies of direct imaging and SP ADE are assumed to be the same, meaning that [ 25 ] ∞ X q =0 H ( q | y ) = τ , (4.5) where τ is the same factor as that for direct imaging. Then N = E [ n ( X )] = ¯ n ( X ) = ν (1) = τ Z dF ( y ) , (4.6) the expected photon number receiv ed in total, is also the same. V . MOMENT ESTIMA TION WITH DIRECT IMA GING Consider the direct-imaging model given by Eq. ( 4.2 ). As- sume that H can be expanded in a T aylor series as H ( x − y ) = ∞ X j =0 ( − 1) j j ! ∂ j H ( x ) ∂ x j y j , (5.1) which leads to f ( x ) = ∞ X j =0 ( − 1) j j ! ∂ j H ( x ) ∂ x j θ j , (5.2) where the unknown parameters are the object moments de- fined by θ j = Z y j dF ( y ) , j ∈ N 0 . (5.3) For the CRB to hold, the parameter space should correspond to the condition that F contains an infinite number of point sources with different positions, as discussed in Appendix B . Appendix C shows a way of reconstructing F from θ via an orthogonal-series expansion, follo wing Ref. [ 51 ]. The score function for each θ j is S j ( x ) = ( − 1) j j ! f ( x ) ∂ j H ( x ) ∂ x j . (5.4) It turns out that the tangent space T for this problem is equal to the whole Hilbert space H under certain technical conditions, as shown in Appendix D . Let the parameter of interest be β = u > θ = ∞ X j =0 u j θ j , (5.5) where u is independent of θ . T o find a candidate influence function, a trick [ 63 ] is to consider the image moments φ gi ven by φ = Z ˜ φ ( x ) d ¯ n ( x ) = ν ( ˜ φ ) , (5.6) where ˜ φ j ( x ) = x j , j ∈ N 0 (5.7) are the monomials. Assuming that all the moments of F and H are finite such that all the moments of f are also finite, φ can be related to the object moments via φ j = Z Z x j H ( x − y ) dF ( y ) dx (5.8) = Z Z ( z + y ) j H ( z ) dF ( y ) dz (5.9) = Z Z j X k =0  j k  z j − k y k H ( z ) dF ( y ) dz (5.10) = ∞ X k =0 C j k θ k , (5.11) where C j k = 1 j ≥ k  j k  Z H ( x ) x j − k dx (5.12) and 1 proposition = ( 1 if proposition is true , 0 otherwise . (5.13) 5 C is a lo wer-triangular matrix, and with C j j = R H ( x ) dx = τ > 0 , C − 1 is well defined and also lo wer-triangular even if the dimension of C is infinite, as shown in Appendix E . The object moments can then be related to the image moments by θ = C − 1 φ, (5.14) and β can be expressed as β = u > θ = u > C − 1 φ = ν ( u > C − 1 ˜ φ ) . (5.15) According to Eq. ( 3.15 ), an influence function is ˜ β ( x ) = u > C − 1 ˜ φ ( x ) = u > ˜ θ ( x ) , (5.16) ˜ θ ( x ) = C − 1 ˜ φ ( x ) . (5.17) Since T = H as shown in Appendix D , the ˜ β given by Eq. ( 5.16 ) belongs to H = T and is efficient according to Eq. ( 3.9 ) as long as it is square-summable. For example, if u contains a finite number of nonzero entries, ˜ β is a polynomial of x and must be square-summable, since all the moments of f are assumed to be finite. The CRB is hence CRB (direct) = ν ( ˜ β 2 ) = u > ν ( ˜ θ ˜ θ > ) u = u > C − 1 ν ( ˜ φ ˜ φ > ) C −> u, (5.18) where C −> = ( C − 1 ) > . This result coincides with that de- riv ed in Ref. [ 25 ] via a more direct but less rigorous method, which is repeated in Appendix F for completeness. An unbiased and efficient estimator ˇ β ( n ) in terms of the observed process n can be constructed from the efficient in- fluence as ˇ β ( n ) = Z ˜ β ( x ) dn ( x ) = u > ˇ θ ( n ) , (5.19) where the object moment estimator is ˇ θ ( n ) = C − 1 ˇ φ ( n ) , ˇ φ ( n ) = Z ˜ φ ( x ) dn ( x ) . (5.20) ˇ β ( n ) is a linear filter of n , so its variance is V ( ˇ β ) = ν ( ˜ β 2 ) , which coincides with the CRB. It is important to note that this estimator does not require any knowledge of the unkno wn pa- rameters, as ˇ φ ( n ) is simply the empirical moments of the ob- served image, and C − 1 is a lower -triangular matrix that de- pends on the moments of the point-spread function H . The estimator still works even if the object happens to consist of a finite number of point sources and θ is on the boundary of the parameter space, although its efficienc y in that case is a more difficult question, as explained in Appendix B . Unlike some of the previous studies on superresolution [ 17 – 20 ], the results here place no restriction on the separations of the point sources and also account for Poisson noise explicitly . VI. CONSTRAINED CRAM ´ ER-RA O BOUND In imaging, the parameters of interest are often the mo- ments with respect to a normalized object distribution with R dF ( y ) = 1 . A simple way of modeling this is to assume that θ 0 = 1 is kno wn. This constraint also makes the model comparable to those in Refs. [ 24 , 26 , 51 ]. Then N = φ 0 = τ θ 0 (6.1) is kno wn as well, implying the constraint γ ( f ) = R f ( x ) dx − N = 0 . The constraint can be differentiated to yield ˙ γ ( f , S j ) = ν ( S j ) = h S j , 1 i = 0 , leading to T ⊥ = span(1) . The new efficient influence, according to Eqs. ( 3.19 ) and ( 3.20 ), should therefore be ˜ β eff = ˜ β − Π( ˜ β |T ⊥ ) = ˜ β − ν ( ˜ β ) ν (1) = ˜ β − β N . (6.2) The constrained CRB is now CRB (direct) θ 0 =1 = ν ( ˜ β 2 eff ) = 1 N h ν 0 ( ˜ β 2 0 ) − β 2 i , (6.3) ˜ β 0 ( x ) = N ˜ β ( x ) = u > ( C /τ ) − 1 ˜ φ ( x ) , (6.4) where ν 0 ( h ) = ν ( h ) /ν (1) is the normalized version of ν . The CRB is necessarily lo wered by the constraint. Other ap- proaches to the constrained CRB yield the same result, as dis- cussed in Appendix G . T o construct a near-ef ficient estimator , suppose that n ( X ) = R dn ( x ) = L > 0 photons have been de- tected. Then dn ( x ) = P L l =1 1 x = X l , and the photon posi- tions { X 1 , X 2 , . . . , X L } are independent and identically dis- tributed according to the probability measure ¯ n/ N . Consider the estimator ˇ β ( n ) = 1 L Z ˜ β 0 ( x ) dn ( x ) = 1 L L X l =1 ˜ β 0 ( X l ) . (6.5) It is straightforward to sho w that E ( ˇ β | n ( X ) = L ) = ν 0 ( ˜ β 0 ) = β , (6.6) V ( ˇ β | n ( X ) = L ) = 1 L h ν 0 ( ˜ β 2 0 ) − β 2 i , (6.7) which is close to the CRB giv en by Eq. ( 6.3 ) if L is close to its expected v alue N . The results are then consistent with standard results in semiparametric estimation concerning the moments of a normalized probability measure [ 21 ]. VII. EVEN-MOMENT ESTIMA TION WITH SP ADE Now consider the SP ADE model giv en by Eqs. ( 4.3 ) and ( 4.4 ) and the Fourier transforms Ψ( k ) = 1 √ 2 π Z ψ ( z ) exp( − ikz ) dz , (7.1) Φ q ( k ) = 1 √ 2 π Z κ ∗ ( q , z ) exp( − ik z ) dz . (7.2) Suppose that Φ = { Φ q ( k ) } is the P AD basis [ 25 , 36 ] given by Φ q ( k ) = √ τ b q ( k )Ψ( k ) , q ∈ N 0 , (7.3) 6 where b = { b q ( k ) : q ∈ N 0 } is the set of orthonormal poly- nomials defined by Z | Ψ( k ) | 2 b q ( k ) b r ( k ) dk = δ q r . (7.4) The polynomials e xist for all q ∈ N 0 as long as the support of | Ψ( k ) | 2 is infinite [ 64 ], and the orthonormality of Φ ensures that the measurement can be implemented by passiv e linear optics [ 23 , 30 , 36 ]. Equation ( 4.4 ) becomes H ( q | y ) = τ     Z | Ψ( k ) | 2 b q ( k ) exp( − ik y ) dk     2 (7.5) = τ       Z | Ψ( k ) | 2 b q ( k ) ∞ X j =0 ( − ik y ) j j ! dk       2 . (7.6) As the b polynomials are derived by applying the Gram- Schmidt procedure to the monomials (1 , k , k 2 , . . . ) > , their basic properties include R | Ψ( k ) | 2 b q ( k ) k r dk = 0 if r < q , R | Ψ( k ) | 2 b q ( k ) k q dk 6 = 0 , and b q ( k ) = ( − 1) q b q ( − k ) if | Ψ( k ) | 2 is ev en, as is often the case in optics. These prop- erties lead to H ( q | y ) = ∞ X j =0 C q j y 2 j , (7.7) where C is an upper-triangular matrix ( C q j = 0 if j < q ) with positiv e diagonal entries ( C q q > 0 ). Equation ( 4.3 ) becomes f ( q ) = ∞ X j =0 C q j θ 2 j , (7.8) which depends on the ev en moments θ 2 j = Z y 2 j dF ( y ) , j ∈ N 0 . (7.9) The score function with respect to each θ 2 j becomes S j ( q ) = 1 f ( q ) ∂ f ( q ) ∂ θ 2 j = C q j f ( q ) . (7.10) Appendix H prov es that T = span( S ) = H . T o find a candidate influence function, suppose that Eq. ( 7.8 ) can be in verted to giv e θ 2 j = ∞ X q =0 ( C − 1 ) j q f ( q ) . (7.11) An influence function for β = u > θ according to Eq. ( 3.15 ) is therefore ˜ β ( q ) = u > ˜ θ ( q ) , ˜ θ 2 j ( q ) = ( C − 1 ) j q . (7.12) Since T = H , this ˜ β belongs to T and is efficient as long as it is square-summable. The CRB is hence CRB (SP ADE) = ν ( ˜ β 2 ) = u > ν ( ˜ θ ˜ θ > ) u = u > C − 1 D C −> u, (7.13) D j k = f ( j ) δ j k . (7.14) A more direct but heuristic way of deriving Eq. ( 7.13 ) is shown in Appendix I . An unbiased and efficient estimator in terms of the observed detector counts n is ˇ β ( n ) = ∞ X q =0 ˜ β ( q ) n ( q ) = u > ∞ X q =0 ˜ θ ( q ) n ( q ) . (7.15) This estimator has a variance V ( ˇ β ) = ν ( ˜ β 2 ) = CRB (SP ADE) , requires no knowledge of any unknown parameter, and still works e ven if the object happens to consist of a finite number of point sources, with no restriction on their separations. If θ 0 = 1 , the constrained CRB can be derived in ways similar to Sec. VI and Appendix G . T o estimate the odd moments of F via SP ADE, variations of the P AD basis are needed [ 24 , 25 ]. The model is much more complicated and a deriv ation of the CRB and the effi- cient estimator is too tedious to work out here, but the upshot is the same: it can be shown that the tangent space for the problem encompasses the whole Hilbert space H , the efficient influence can be retrieved from the relation β = ν ( ˜ β ) , and an unbiased and efficient estimator is ˇ β ( n ) = R ˜ β ( x ) dn ( x ) . A. Gaussian point-spread function More explicit results can be obtained and the assumptions can be checked more carefully by assuming the Gaussian point-spread function ψ ( z ) = 1 (2 π ) 1 / 4 exp  − z 2 4  . (7.16) The P AD basis becomes the Hermite-Gaussian basis, and it can be shown that [ 23 , 24 , 55 ] H ( q | y ) = τ exp  − y 2 4  ( y / 2) 2 q q ! . (7.17) The C matrix in Eq. ( 7.7 ) can be determined by expanding exp( − y 2 / 4) = P ∞ j =0 ( − y 2 / 4) j /j ! , gi ving C q j = 1 j ≥ q τ ( − 1) j − q 4 j q !( j − q )! . (7.18) It is not difficult to check that the matrix in verse of C is ˜ θ 2 j ( q ) = ( C − 1 ) j q = 1 q ≥ j 4 j q ! τ ( q − j )! , (7.19) which is a degree- j polynomial of q . P ∞ q =0 ˜ θ 2 j ( q ) f ( q ) is the j th factorial moment of f and indeed equal to θ 2 j , since H ( q | y ) is Poisson and its factorial moment is P ∞ q =0 ˜ θ 2 j ( q ) H ( q | y ) = y 2 j [ 65 ]. In general, each degree- j moment of H ( q | y ) is a degree- j polynomial of y 2 , so each degree- j moment of f ( q ) is a linear combination of the mo- ments of F up to degree 2 j . All the moments of f are there- fore finite as long as all the moments of F are finite. If u has a finite number of nonzero entries, the influence function giv en by Eqs. ( 7.12 ) is a polynomial of q , so ν ( ˜ β 2 ) < ∞ , and ˜ β ∈ H is ensured. 7 B. Bandlimited point-spread function Another important example is the bandlimited point-spread function giv en by Ψ( k ) = 1 | k | < 1 √ 2 . (7.20) b is then the well known set of Legendre polynomials [ 66 ]. Appendix J shows the detailed calculations; here I list the re- sults only . Equation ( 7.5 ) becomes H ( q | y ) = τ (2 q + 1) j 2 q ( y ) , (7.21) where j q ( y ) is the spherical Bessel function of the first kind [ 66 , Eq. (10.47.3)]. An influence function for estimating θ 2 j with θ 2 j = ν ( ˜ θ 2 j ) is ˜ θ 2 j ( q ) = 1 q ≥ j (2 j + 1)!!(2 j − 1)!! τ  q + j 2 j  , (7.22) where !! denotes the double factorial [ 66 ]. ˜ θ 2 j ( q ) is a degree- 2 j polynomial of q , so ˜ β ( q ) is also a polynomial of q if u contains a finite number of nonzero entries. As long as all the moments of F are finite, all the moments of f can also be shown to be finite, and ν ( ˜ β 2 ) < ∞ is ensured. Notice that the direct-imaging theory in Sec. V breaks down for this bandlimited point-spread function, as the second and higher ev en moments of H ( x ) = τ | ψ ( x ) | 2 = ( τ /π ) sinc 2 ( x ) are all infinite. The CRB in that case remains an open prob- lem, although it is possible to apodize the point-spread func- tion optically such that all its moments become finite and the semiparametric estimator giv en by Eq. ( 5.19 ) has a finite vari- ance. For example, if Ψ( k ) ∝ 1 | k | < 1 exp  − 1 k 2 − 1  , (7.23) then Ψ( k ) is infinitely differentiable despite the hard band- width limit [ 67 ] and all the moments of | ψ ( x ) | 2 are finite [ 25 ]. VIII. COMP ARISON OF IMA GING METHODS The adv antage of SP ADE o ver direct imaging occurs in the subdiffraction regime, where the width ∆ of the object dis- tribution F with respect to the origin is much smaller than the width of the point-spread function ψ [ 24 – 26 , 51 ]. As the width of ψ is normalized as 1 , the regime is defined as ∆  1 , (8.1) and the object moments scale as θ j = θ 0 O (∆ j ) . (8.2) W ith the attainable CRBs giv en by Eqs. ( 5.18 ) and ( 7.13 ) at hand, I can now compare direct imaging and SP ADE on the same semiparametric footing. Consider the estimation of a specific moment θ k with u j = δ j k . (8.3) For direct imaging in the subdiffraction re gime, the image be- comes close to the point-spread function, viz., f ( x ) ≈ θ 0 H ( x ) = N | ψ ( x ) | 2 , (8.4) where N , the expected photon number recei ved in total, is giv en by Eq. ( 4.6 ). W ith C j k = τ O (1) and ν ( ˜ φ ˜ φ > ) = N O (1) , the CRB becomes CRB (direct) = θ 2 0 N O (1) . (8.5) For SP ADE on the other hand, notice that the C and C − 1 matrices are upper-triangular , meaning that f ( q ) = N O (∆ 2 q ) , (8.6) and the CRB for estimating θ k , where k is e ven, becomes CRB (SP ADE) = θ 2 0 N O (∆ k ) , (8.7) which is much lower than Eq. ( 8.5 ) when ∆  1 and k ≥ 2 . This is consistent with earlier approximate results [ 24 , 25 ]. An intuitive explanation of the enhancement, as well as a discussion of the limitations of SP ADE, can be found in Ref. [ 51 ]. The constrained CRB with θ 0 = 1 becomes [ O (∆ k ) − θ 2 k ] / N = O (∆ k ) / N , which is on the same order of magnitude as the fundamental quantum limit [ 26 ]. More exact and pleasing results can be obtained if ψ is Gaussian and given by Eq. ( 7.16 ). Consider for example the estimation of the second moment θ 2 . For direct imaging, it can be shown that CRB (direct) = 1 τ (2 θ 0 + 4 θ 2 + θ 4 ) = θ 2 0 N O (1) . (8.8) For SP ADE on the other hand, CRB (SP ADE) = 1 τ (4 θ 2 + θ 4 ) = θ 2 0 N O (∆ 2 ) , (8.9) which not only beats direct imaging by a significant amount in the subdiffraction regime but is in fact superior for all parame- ter v alues. T o further illustrate the dif ference between the tw o methods, suppose that the object happens to be a flat top gi ven by dF ( y ) = θ 0 ∆ 1 | y | < ∆ / 2 dy . (8.10) Do note that the semiparametric CRBs do not assume the knowledge of this object shape, which is specified here only for the purpose of plotting examples of the CRBs. With θ 2 = θ 0 ∆ 2 / 12 and θ 4 = θ 0 ∆ 4 / 80 , Fig. 3 plots Eqs. ( 8.8 ) and ( 8.9 ) against ∆ in log-log scale. The g ap between the two curves in the ∆  1 re gime is striking. 8 FIG. 3. The semiparametric CRBs for the second moment θ 2 giv en by Eqs. ( 8.8 ) and ( 8.9 ) versus the object size ∆ in log-log scale, if the point-spread function is Gaussian and the object happens to be a flat top. Both the CRBs and ∆ are normalized such that they are dimensionless. W ith the constraint θ 0 = 1 , the CRBs become CRB (direct) θ 0 =1 = 1 N  2 + 4 θ 2 + θ 4 − θ 2 2  = O (1) N , (8.11) CRB (SP ADE) θ 0 =1 = 1 N  4 θ 2 + θ 4 − θ 2 2  = O (∆ 2 ) N . (8.12) It is note worthy that Eq. ( 8.12 ) is exactly equal to the quan- tum limit giv en by [ 51 , Eq. (E15)], meaning that SP ADE is exactly quantum-optimal—at all parameter values—for esti- mating the second moment. This is consistent with previous results concerning the estimation of two-point separation [ 23 ] and object size [ 24 , 28 ], but note that the pre vious results as- sume that the object shape is kno wn, whereas the CRBs and the estimators here assume the opposite. IX. CONCLUSION The semiparametric theory set forth solves an important and difficult problem in incoherent optical imaging: the ev alu- ation of the CRB and the ef ficient estimation of object parame- ters when little prior information about the object is available. The theory giv es exact and achiev able semiparametric CRBs for both direct imaging and SP ADE, establishing the superi- ority and versatility of SP ADE beyond the special parametric scenarios considered by previous studies. Despite the elegant results, the theory has a fe w shortcom- ings. On the mathematical side, the conditions for the theory to hold seem difficult to check in the case of direct imaging with a non-Gaussian point-spread function, especially when the point-spread function has infinite moments. It is an open question whether this is merely a technicality or a hint at a whole new regime of statistics. On the practical side, the theory may be accused of assuming ideal conditions for both measurements, such as infinitesimal pix els for direct imaging, the av ailability of infinitely many modes for SP ADE, perfect specification and knowledge of the optical systems, and the absence of excess noise. Reality is necessarily uglier , but the results here remain rele vant by serving as fundamental lim- its (via the data-processing inequality [ 51 , 68 ]) and offering insights into the essential physics. The theoretical and exper - imental progress on SP ADE and related methods so far [ 23 – 58 , 69 – 71 ] has provided encouragement that the theory has realistic potential, and the general results here should moti- vate further research into the wider applications of quantum- inspired imaging methods. An interesting future direction is to generalize the semi- parametric formalism for quantum estimation [ 72 , 73 ]. By treating the symmetric logarithmic deri vati ves of the quantum state ρ as the scores in the L 2 h ( ρ ) space proposed by Holev o [ 73 ] and adopting a geometric picture [ 74 ], a quantum gen- eralization of the semiparametric CRB can be envisioned, b ut whether it can solv e any important problem, such as the quan- tum limit to incoherent imaging [ 26 , 27 ], remains to be seen. A CKNO WLEDGMENTS This work is supported by the Singapore National Research Foundation under Project No. QEP-P7. Appendix A: Pr oof of the semiparametric CRB for Poisson processes Define the inner product between two random variables ˇ r 1 and ˇ r 2 as ( ˇ r 1 , ˇ r 2 ) = E ( ˇ r 1 ˇ r 2 ) , (A1) and the norm as ||| ˇ r ||| = p ( ˇ r , ˇ r ) = p E ( ˇ r 2 ) . (A2) Let the Hilbert space of zero-mean random variables be ˇ R =  ˇ r : E ( ˇ r ) = 0 , E ( ˇ r 2 ) < ∞  , (A3) and define ˇ T = span( ˇ S ) ⊆ ˇ R , (A4) where ˇ S is defined by Eq. ( 2.4 ). Let ˇ δ ∈ ˇ R be any random variable that satisfies E ( ˇ δ ˇ S ) = u. (A5) The semiparametric CRB is [ 21 , 22 ] E ( ˇ δ 2 ) ≥ CRB = E ( ˇ δ 2 eff ) , (A6) ˇ δ eff = Π( ˇ δ | ˇ T ) = arg min ˇ h ∈ ˇ T ||| ˇ δ − ˇ h ||| . (A7) The proof can be done via a Pythagorean theorem [ 22 ] without recourse to the Cauchy-Schwartz inequality or the existence 9 of J − 1 . ˇ S j is called a score and ˇ δ an influence in statistics [ 21 , 22 ]; this paper uses the same terminology for S and ˜ β in light of their resemblance to the statistical quantities. The resemblance can be turned into a precise correspon- dence for a Poisson random measure by considering the sub- space ˇ H ⊂ ˇ R of random v ariables that are linear with respect to n . Any element ˇ h ∈ ˇ H can be expressed as ˇ h = U h = Z h ( x ) [ dn ( x ) − d ¯ n ( x )] , (A8) where U is a surjectiv e linear map from H to ˇ H . Since ( U h 1 , U h 2 ) = h h 1 , h 2 i ∀ h 1 , h 2 ∈ H (A9) by virtue of Eq. ( 2.3 ), ˇ H is isomorphic to H and U is unitary [ 61 ], and since ˇ T ⊆ ˇ H and ˇ S = U S , ˇ T is isomorphic to T . Picking a ˇ δ ∈ ˇ H with ˇ δ = U ˜ β = Z ˜ β ( x ) [ dn ( x ) − d ¯ n ( x )] (A10) leads to E ( ˇ δ ˇ S ) = ν ( ˜ β S ) = u, ˇ δ eff = U ˜ β eff , (A11) where ˜ β eff is giv en by Eq. ( 3.9 ) because of Eq. ( A7 ) and the isomorphisms. The CRB becomes CRB = E ( ˇ δ 2 eff ) = ν ( ˜ β 2 eff ) , (A12) which is Eq. ( 3.8 ). Appendix B: The moment parameter space Define an s × s Hankel matrix with respect to a real-number sequence θ = ( θ 0 , θ 1 , . . . ) > as M ( s ) j k ( θ ) = θ j + k , j, k ∈ { 0 , 1 , . . . , s − 1 } . (B1) If θ is a moment sequence that arises from a nonnegati ve mea- sure F , w > M ( s ) w = Z   s − 1 X j =0 w j y j   2 dF ( y ) (B2) is nonne gati ve for any real vector w , and all Hankel matrices are positiv e-semidefinite, viz., M ( s ) ≥ 0 ∀ s ∈ N 1 . (B3) Con versely , any sequence with Hankel matrices that obey Eq. ( B3 ) can be e xpressed in the form of Eq. ( 5.3 ) with a non- negati ve F by virtue of Hambur ger’ s theorem [ 75 ]. For the CRB to hold for a p -dimensional θ , the parameter space Θ should be an open subset of R p [ 68 , 76 ]. Intuitively , the requirement makes sense because all the parameters in θ are unknown and θ should be allowed to vary in a neighbor- hood, otherwise the problem would be overparametrized. If Θ is not an open subset, the parameter space would be con- strained and the CRB may be modified [ 76 ]. When all the moments are unknown parameters, consider the set Θ = n θ : M ( s ) ( θ ) > 0 ∀ s ∈ N 1 o . (B4) Each θ ∈ Θ corresponds to a measure with infinite support r = ∞ [ 75 ]. The proof can be done by observing that the polynomial in Eq. ( B2 ) has at most s − 1 zeros and the integral is strictly positiv e for any w 6 = 0 if and only if r ≥ s , and therefore the constraint for Θ is satisfied if and only if r = ∞ . For r < ∞ , F can be expressed in terms of its support { y l : 0 ≤ l ≤ r − 1 , y l < y l +1 } as dF ( y ) = r − 1 X l =0 F l 1 y = y l , dF ( y ) dy = r − 1 X l =0 F l δ ( y − y l ) . (B5) In the context of optics, r is the minimum number of point sources that can describe the object distrib ution. The assump- tion of Eq. ( B4 ) as the parameter space is consistent with the infinite-support assumption for semiparametric estimation with mixture models [ 21 , Sec. 6.5], and it also makes intuitiv e sense, at least as a necessary condition—with r point sources, there are only 2 r unknown parameters, and the problem would be ov erparametrized if all the moments are assumed to be un- known. Further inequality constraints on θ may be needed to ensure the con ver gence of the T aylor series in Eqs. ( 5.2 ) and ( 7.8 ), although they should not affect the CRB as long as θ obeys and stays a way from them [ 76 ]. The boundary of Θ corresponds to measures with finite sup- port r < ∞ . If s ≤ r , then M ( s ) > 0 and M ( s ) is full-rank (rank = s ), but if s > r , I can write M ( s ) = V > diag( F ) V , (B6) V j k = y k l , diag( F ) j k = 1 0 ≤ j ≤ r − 1 F j δ j k . (B7) V is the V andermonde matrix and in vertible since { y l } are assumed to be distinct [ 77 ]. With M ( s ) ≥ 0 and diag( F ) ≥ 0 , Sylvester’ s law of inertia [ 77 ] implies that the rank of M ( s ) is the same as the rank of diag( F ) , which is r . In other words, the rank of M ( s ) is min( r , s ) , and any finite r means that M ( s ) is rank-deficient and does not satisfy the strict M ( s ) > 0 for all s > r . Whether the CRB still holds for θ on the boundary is a difficult question, considering that the parameter space here is infinite-dimensional and it is not obvious how existing finite-dimensional results regarding the CRB on a boundary [ 76 ] can be applied. Appendix C: Series expansion of the object distrib ution Assume that the object measure F can be expressed as the orthogonal series dF ( y ) = ∞ X j =0 ξ j g j ( y ) dF (0) ( y ) (C1) 10 with respect to a reference measure F (0) , where { g j = P ∞ k =0 G j k y k : j ∈ N 0 } are the orthogonal polynomials de- fined by g j ( y ) = ∞ X k =0 G j k y k , Z g j ( y ) g k ( y ) dF (0) ( y ) = δ j k , (C2) and G is a lower -triangular matrix with nonzero diagonal en- tries that can be obtained by the Gram-Schmidt procedure. Thus each “Fourier” coefficient ξ j can be expressed in terms of the moment parameters as ξ j = Z g j ( y ) dF ( y ) = ∞ X k =0 G j k θ k , (C3) which can be written as ξ = Gθ , θ = G − 1 ξ . (C4) Hence each θ corresponds to a set of coefficients ξ that can be used to represent F via Eq. ( C1 ). It is straightforward to transform the CRBs and the efficient estimators deriv ed in this paper for θ to those for ξ via Eqs. ( C4 ). Appendix D: T angent space f or the direct-imaging model Consider the tangent space T given by Eq. ( 3.10 ) and the score functions given by Eq. ( 5.4 ) for direct imaging. First note that S ⊂ H , as the Fisher information J j j = h S j , S j i = ν ( S 2 j ) is assumed to be finite for all j . Recent calculations in quantum estimation theory suggest that this assumption is rea- sonable for any measurement [ 26 ]. T o prove T = span( S ) = H , the standard method is to show that the only element in H orthogonal to S is 0 [ 61 ], that is, h h, S j i = 0 ∀ j ∈ N 0 (D1) implies h = 0 (almost ev erywhere with respect to ¯ n ). Here I list a few approaches with v arious lev els of rigor . The first approach is to consider the set of orthogonal poly- nomials a = n a j ( x ) = A ˜ φ ( x ) : j ∈ N 0 , h a j , a k i = δ j k o , (D2) where A is a lower-triangular matrix with nonzero diagonal entries and can be determined by applying the Gram-Schmidt procedure to the monomials ˜ φ ( x ) [ 64 ]. Under rather general conditions on f , the polynomials form an orthonormal basis of H [ 64 ], viz., H = span( a ) , (D3) and I can write Eq. ( D1 ) as h h, S j i = ∞ X k =0 h h, a k i h a k , S j i = 0 ∀ j ∈ N 0 , (D4) or , more compactly , B > w = 0 , w k = h h, a k i , B kj = h a k , S j i . (D5) If the only solution to Eq. ( D5 ) is w = 0 , then the only solu- tion to Eq. ( D4 ) is h = 0 . This is equiv alent to the condition that B > is injectiv e. Integration by parts yields B kj = ( − 1) j j ! Z a k ( x ) ∂ j H ( x ) ∂ x j dx = ∞ X l =0 A kl C lj , (D6) where C is the same as Eq. ( 5.12 ). Since both A and C are lower -triangular with nonzero diagonal entries, B = AC is also lower-triangular with nonzero diagonal entries, and B > has a well defined matrix inv erse ( B > ) − 1 = ( B − 1 ) > = A −> C −> in the sense that B > ( B > ) − 1 = ( B > ) − 1 B > = I , (D7) where I is the identity matrix; see Appendix E for details. If the matrices were finite-dimensional, the existence of a matrix in verse would imply ( B > ) − 1 ( B > w ) = [( B > ) − 1 B > ] w = w , (D8) and the only solution to Eq. ( D5 ) w ould be w = 0 . This proof is not entirely satisfactory howe ver , as Eq. ( D8 ) assumes that the product of the infinite-dimensional matrices is associativ e. Associativity assumes that the order of the sums inv olved in the matrix product can be interchanged, but it cannot be guar - anteed for infinite-dimensional matrices. In other words, the existence of a matrix in verse for B > may not imply that B > is injectiv e. A more rigorous approach is to define χ y ( x ) = ∞ X j =0 y j S j ( x ) , y ∈ Y ⊂ R , (D9) and notice that Eq. ( D1 ) implies h h, χ y i = Z h ( x ) ∞ X j =0 y j ( − 1) j j ! ∂ j H ( x ) ∂ x j dx (D10) = Z h ( x ) H ( x − y ) dx = 0 ∀ y ∈ Y . (D11) The unique solution to Eq. ( D11 ) is h = 0 if the family { H ( x − y ) : y ∈ Y } satisfies a property called complete- ness in statistics [ 5 ]. For example, if H is Gaussian, { H } is a full-rank exponential family for any open subset Y ⊂ R and therefore complete [ 5 ]. An even more rigorous formulation of this approach [ 21 ] is to treat h h, χ y i as an operator that maps h ∈ H to a function of y in another Hilbert space, and then sho w that the null space of the operator consists of only h = 0 . The proof again boils down to the requirement that { H } should be complete; see Ref. [ 21 , Sec. 6.5]. 11 Appendix E: In verse of an infinite-dimensional triangular matrix Let C be an infinite-dimensional matrix with entries in- dex ed by ( j , k ) ∈ N 2 0 . Define its formal matrix in verse C − 1 as another infinite-dimensional matrix that satisfies ∞ X l =0 C j l ( C − 1 ) lk = δ j k . (E1) If C is lo wer-triangular with nonzero diagonal entries, viz., C j k = 0 if k > j, C j j 6 = 0 , (E2) then C − 1 can be found by a recursi ve relation as follows. De- fine C ( s ) as the s × s upper-left submatrix of C . Write C ( s +1) and ( C − 1 ) ( s +1) as the partitions C ( s +1) =  C ( s ) 0 c > C ss  , (E3) ( C − 1 ) ( s +1) =  ( C − 1 ) ( s ) 0 d > ( C − 1 ) ss  . (E4) Giv en ( C − 1 ) ( s ) = ( C ( s ) ) − 1 , d > = − c > ( C ( s ) ) − 1 C ss , ( C − 1 ) ss = 1 C ss , (E5) and the recursion starts from ( C − 1 ) (1) = ( C (1) ) − 1 with one element ( C − 1 ) 00 = 1 /C 00 . The matrix in verse of an infinite- dimensional upper-triangular matrix can be defined in a simi- lar way . Although the product of infinite-dimensional matrices may not be associati ve, it can still be proved by induction that D ( C u ) = ( D C ) u for any vector u if D and C are lo wer- triangular , because D ( C u ) = ∞ X k =0 D j k ∞ X l =0 C kl u l = j X k =0 D j k k X l =0 C kl u l (E6) in volves finite sums only . Thus it is safe to assume that C − 1 ( C u ) = ( C − 1 C ) u = u and C ( C − 1 u ) = ( C C − 1 ) u = u if C is lo wer-triangular with nonzero diagonal entries. Appendix F: An alter native derivation of the Cram ´ er -Rao bound for dir ect imaging Consider the problem described in Sec. V . Since the poly- nomials given by Eq. ( D2 ) are an orthonormal basis, the infor - mation matrix for the moment parameters can be e xpressed as J j k = h S j , S k i = ∞ X l =0 h S j , a l ih a l , S k i , (F1) meaning that J = B > B , where B = AC is giv en by Eq. ( D6 ). Ignoring the complications of multiplying infinite- dimensional matrices, the CRB becomes J − 1 = B − 1 B −> = C − 1 A − 1 A −> C −> . (F2) T o ev aluate A − 1 A −> , notice that the orthonormality of a can be written as h a j , a k i = ∞ X l =0 ∞ X m =0 A j l D ˜ φ l , ˜ φ m E A km = δ j k , (F3) where ˜ φ is the monomials gi ven by Eq. ( 5.7 ). In other words, Aν ( ˜ φ ˜ φ > ) A > = I , A − 1 A −> = ν ( ˜ φ ˜ φ > ) , (F4) giving J − 1 = C − 1 ν ( ˜ φ ˜ φ > ) C −> . (F5) This leads to Eq. ( 5.18 ) for the parameter β = u > θ . Appendix G: Alter native approaches to the constrained Cram ´ er -Rao bound One way of deri ving the constrained CRB if θ 0 is kno wn is to consider the information matrix ˜ J with respect to the pa- rameters ϑ = ( θ 1 , θ 2 , . . . ) > without θ 0 . Then θ = ( θ 0 , ϑ > ) > , and ˜ J can be written as the submatrix of J , or J =  J 00 j > j ˜ J  , (G1) where j is a column v ector . Ignore the complications of deal- ing with infinite-dimensional matrices and partition J − 1 sim- ilarly as J − 1 =  E 00 e > e ˜ E  . (G2) Then it is straightforward to sho w that ˜ J − 1 = ˜ E − ee > E 00 . (G3) Let ˜ ϑ = ( ˜ θ 1 , ˜ θ 2 , . . . ) > , and observe that ˜ θ 0 = 1 /C 00 from Eqs. ( 5.17 ), ( 5.12 ), and ( 5.6 ). Then Eq. ( 5.18 ) implies that ˜ E = ν ( ˜ ϑ ˜ ϑ > ) , (G4) e = ν ( ˜ ϑ ˜ θ 0 ) = ν ( ˜ ϑ ) C 00 = ϑ C 00 , (G5) E 00 = ν ( ˜ θ 0 ˜ θ 0 ) = ν (1) C 2 00 = φ 0 C 2 00 . (G6) Hence ˜ J − 1 = ν ( ˜ ϑ ˜ ϑ > ) − ϑϑ > φ 0 , (G7) which implies Eq. ( 6.3 ) if the parameter of interest is defined as β = u > θ with u 0 = 0 . Y et another way of deriving the constrained CRB can be found in Ref. [ 76 ], which can be shown to lead to the same result here. 12 Appendix H: T angent space f or the SP ADE model The proof is similar to the first approach in Appendix D . Consider H = span( a ) in terms of an obvious orthonormal basis a = n a j ( q ) = δ j q / p f ( j ) : j ∈ N 0 o . (H1) Any h ∈ H orthogonal to the S giv en by Eq. ( 7.10 ) obeys h h, S j i = ∞ X k =0 h h, a k i h a k , S j i = 0 ∀ j ∈ N 0 , (H2) which can be written as B > w = 0 , w k = h h, a k i , (H3) B j k = h a j , S k i = C j k p f ( j ) . (H4) As C is upper-triangular with nonzero diagonal entries and f > 0 is assumed, B > is lower -triangular with nonzero di- agonal entries, and induction can be used to prov e that the only solution to B > w = 0 is w = 0 , or equiv alently h = 0 . Hence T = span( S ) = H . The proof is easier than the one in Appendix D because B > here is lower -triangular rather than upper-triangular . An alternati ve proof, similar to the second approach in Ap- pendix D and Ref. [ 21 , Sec. 6.5] b ut less fruitful in this case, is to define χ y ( x ) = ∞ X j =0 y 2 j S j ( x ) , y ∈ Y ⊂ R , (H5) consider h h, χ y i = ∞ X q =0 h ( q ) H ( q | y ) = 0 , (H6) and use the completeness of { H ( q | y ) : y ∈ Y } to prov e the unique solution h = 0 . If H is Poisson, for example, then { H } is a full-rank e xponential family and therefore complete [ 5 ]. Appendix I: An alter native deri vation of the Cram ´ er -Rao bound for SP ADE Consider the problem described in Sec. VII . W ith the or- thonormal basis giv en by Eq. ( H1 ) and the B matrix giv en by Eq. ( H4 ), the information matrix with respect to the moment parameters can again be e xpressed as J = B > B according to Eq. ( F1 ). With Eq. ( H4 ), B − 1 becomes ( B − 1 ) j q = ( C − 1 ) j q p f ( q ) . (I1) Ignoring the complications of multiplying infinite- dimensional matrices, the CRB J − 1 = B − 1 B −> is ( J − 1 ) j k = ∞ X q =0 ( C − 1 ) j q f ( q )( C − 1 ) kq = C − 1 D C −> , (I2) where D is giv en by Eq. ( 7.14 ), and the CRB for β = u > θ coincides with Eq. ( 7.13 ). Appendix J: Calculations concer ning SP ADE for a bandlimited point-spread function Consider the point-spread function gi ven by Eq. ( 7.20 ). The standard Legendre polynomials are defined in terms of 1 2 Z 1 − 1 P q ( k ) P p ( k ) dk = 1 2 q + 1 δ q p , (J1) such that the orthonormal version is b q ( k ) = p 2 q + 1 P q ( k ) . (J2) The Fourier transform of the polynomial is [ 66 , Eq. (18.17.19)] 1 2 Z 1 − 1 b q ( k ) exp( ik y ) dk = p 2 q + 1 j q ( y ) , (J3) where j q ( y ) is the spherical Bessel function of the first kind [ 66 , Eq. (10.47.3)]. Then Eq. ( 7.21 ) follows from Eq. ( 7.5 ) and ( J3 ). Let ˜ H ( q | y ) = H ( q | y ) τ = (2 q + 1) j 2 q ( y ) . (J4) From Ref. [ 66 , Eq. (10.60.2)], one can deriv e the useful for- mula ∞ X q =0 ˜ H ( q | y ) P q ( k ) = sinc w = ( (sin w ) /w , w 6 = 0 , 1 , w = 0 , (J5) w = y √ 2 − 2 k . (J6) For example, since P q (1) = 1 , one can check that P ∞ q =0 ˜ H ( q | y ) = 1 in accordance with Eq. ( 4.5 ). Using the facts sinc w = 1 2 Z 1 − 1 exp( iw z ) dz = ∞ X l =0 ( − 1) l w 2 l (2 l + 1)! , (J7) dw dk = − y w , P ( j ) q (1) = d j P q ( k ) dk j     k =1 , (J8) it can also be shown that ∞ X q =0 ˜ H ( q | y ) P ( j ) q (1) = d j sinc w dk j     k =1 = y 2 j (2 j + 1)!! , (J9) which leads to ∞ X q =0 f ( q ) P ( j ) q (1) = τ θ 2 j (2 j + 1)!! . (J10) An influence function for estimating θ 2 j is hence ˜ θ 2 j ( q ) = (2 j + 1)!! τ P ( j ) q (1) , (J11) 13 which obeys θ 2 j = ν ( ˜ θ 2 j ) . T o deriv e an explicit ex- pression for P ( j ) q (1) , consider the Rodrigues formula [ 66 , Eq. (14.7.13)] P q ( k ) = 1 2 q q ! d q dk q ( k 2 − 1) q , (J12) which leads to P q ( k ) = q X l =0  q l   q + l l   k − 1 2  l , (J13) P ( j ) q (1) = 1 q ≥ j (2 j − 1)!!  q + j 2 j  , (J14) and Eq. ( 7.22 ) results. T o bound the moments of ˜ H and f , consider a lower bound on the binomial coefficient for j ≥ 1 given by [ 78 , pp. 1186]  q + j 2 j  ≥ ( q + j ) 2 j (2 j ) 2 j ≥ q 2 j (2 j ) 2 j , (J15) such that each ev en moment of ˜ H can be bounded as ∞ X q =0 ˜ H ( q | y ) q 2 j = j − 1 X q =0 ˜ H ( q | y ) q 2 j + ∞ X q = j ˜ H ( q | y ) q 2 j (J16) ≤ ( j − 1) 2 j + (2 j ) 2 j (2 j − 1)!! ∞ X q =0 ˜ H ( q | y ) P ( j ) q (1) (J17) = ( j − 1) 2 j + (2 j ) 2 j y 2 j (2 j − 1)!!(2 j + 1)!! . (J18) This means that each ev en moment of f ( q ) is bounded as ν ( q 2 j ) ≤ τ  ( j − 1) 2 j θ 0 + (2 j ) 2 j θ 2 j (2 j − 1)!!(2 j + 1)!!  . (J19) W ith ν ( q 0 ) = ν (1) = τ θ 0 , ν ( q 2 j ) < ∞ for all j ∈ N 0 as long as θ 0 and θ 2 j are finite. Odd moments can be bounded via the Cauchy-Schwartz inequality [ ν ( q j )] 2 ≤ ν (1) ν ( q 2 j ) . Hence all the moments of f are finite as long as all the moments of F are finite. [1] Max Born and Emil W olf, Principles of Optics: Electroma g- netic Theory of Pr opagation, Interference and Diffraction of Light (Cambridge University Press, Cambridge, 1999). [2] Joseph W . Goodman, Intr oduction to F ourier Optics (McGraw- Hill, New Y ork, 2004). [3] Leonard Mandel and Emil W olf, Optical Coher ence and Quan- tum Optics (Cambridge Univ ersity Press, Cambridge, 1995). [4] Joseph W . Goodman, Statistical Optics (Wile y , New Y ork, 1985). [5] E. L. Lehmann and George Casella, Theory of P oint Estimation (Springer , New Y ork, 1998). [6] Edward J. Farrell, “Information Content of Photoelectric Star Images, ” Journal of the Optical Society of America 56 , 578– 587 (1966) . [7] Ming-Jer Tsai and K eh-Ping Dunn, P erformance Limitations on P arameter Estimation of Closely Spaced Optical T arg ets Using Shot-Noise Detector Model , T ech. Rep. AD A073462 (Lincoln Laboratory , MIT , 1979). [8] Jonas Zmuidzinas, “Cram ´ er–Rao sensitivity limits for astro- nomical instruments: implications for interferometer design, ” Journal of the Optical Society of America A 20 , 218–233 (2003) . [9] Eric D. Feigelson and G. Jogesh Babu, Modern Statistical Methods for Astr onomy (Cambridge University Press, Cam- bridge, 2012). [10] Sripad Ram, E. Sally W ard, and Raimund J. Ober, “Beyond Rayleigh’ s criterion: A resolution measure with application to single-molecule microscopy , ” Proceedings of the National Academy of Sciences of the United States of America 103 , 4457–4462 (2006) . [11] Alex Small and Shane Stahlheber , “Fluorophore localization al- gorithms for super-resolution microscop y , ” Nature Methods 11 , 267–279 (2014) . [12] Hendrik Deschout, Francesca Cella Zanacchi, Michael Mlodzianoski, Alberto Diaspro, Joerg Be wersdorf, Samuel T . Hess, and Ke vin Braeckmans, “Precisely and accurately local- izing single emitters in fluorescence microscop y , ” Nature Meth- ods 11 , 253–266 (2014) . [13] Jerry Chao, E. Sally W ard, and Raimund J. Ober , “Fisher in- formation theory for parameter estimation in single molecule microscopy: tutorial, ” Journal of the Optical Society of Amer- ica A 33 , B36 (2016) . [14] Alex v on Diezmann, Y oav Shechtman, and W . E. Mo- erner , “Three-Dimensional Localization of Single Molecules for Super-Resolution Imaging and Single-Particle Tracking, ” Chemical Revie ws 117 , 7244–7275 (2017) . [15] E. Bettens, D. V an Dyck, A. J. den Dekker , J. Sijbers, and A. van den Bos, “Model-based two-object resolution from ob- servations having counting statistics, ” Ultramicroscopy 77 , 37– 48 (1999) . [16] S. V an Aert, A. J. den Dekker , D. V an Dyck, and A. v an den Bos, “High-resolution electron microscopy and electron tomog- raphy: resolution versus precision, ” Journal of Structural Biol- ogy 138 , 21–33 (2002) . [17] Geoffrey de V illiers and E. Roy Pike, The Limits of Resolution (CRC Press, Boca Raton, 2016). [18] Emmanuel J. Cand ` es and Carlos Fernandez-Granda, “Super- Resolution from Noisy Data, ” Journal of Fourier Analysis and Applications 19 , 1229–1254 (2013) . [19] Emmanuel J. Cand ` es and Carlos Fernandez-Granda, “T owards a Mathematical Theory of Super-resolution, ” Communications on Pure and Applied Mathematics 67 , 906–956 (2014) . [20] Geoffrey Schiebinger , Elina Robe va, and Benjamin Recht, “Su- perresolution without separation, ” Information and Inference: A Journal of the IMA , iax006 (2017) . [21] Peter J. Bickel, Chris A. J. Klaassen, Y a’acov Ritov , and 14 John A. W ellner, Efficient and Adaptive Estimation for Semi- parametric Models (Springer , New Y ork, 1993). [22] Anastasios A. Tsiatis, Semiparametric Theory and Missing Data (Springer , New Y ork, 2006). [23] Mankei Tsang, Ranjith Nair , and Xiao-Ming Lu, “Quan- tum theory of superresolution for two incoherent optical point sources, ” Physical Re view X 6 , 031033 (2016) . [24] Mankei Tsang, “Subdiffraction incoherent optical imaging via spatial-mode demultiplexing, ” New Journal of Physics 19 , 023054 (2017) . [25] Mankei Tsang, “Subdiffraction incoherent optical imaging via spatial-mode demultiplexing: Semiclassical treatment, ” Physi- cal Revie w A 97 , 023830 (2018) . [26] Mankei Tsang, “Quantum limit to subdiffraction incoherent op- tical imaging, ” Physical Re view A 99 , 012305 (2019) . [27] Sisi Zhou and Liang Jiang, “Modern description of Rayleigh’ s criterion, ” Physical Re view A 99 , 013808 (2019) . [28] Zachary Dutton, Ronan K erviche, Amit Ashok, and Saikat Guha, “ Attaining the quantum limit of superresolution in imag- ing an object’ s length via predetection spatial-mode sorting, ” Physical Re view A 99 , 033847 (2019) . [29] Ranjith Nair and Mankei Tsang, “Interferometric superlocal- ization of tw o incoherent optical point sources, ” Optics Express 24 , 3684–3701 (2016) . [30] Mankei Tsang, Ranjith Nair, and Xiao-Ming Lu, “Quantum information for semiclassical optics, ” in Pr oc. SPIE , Quantum and Nonlinear Optics IV , V ol. 10029 (SPIE, Bellingham, W A, 2016) p. 1002903. [31] Ranjith Nair and Mankei Tsang, “Far -Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit, ” Physical Re view Letters 117 , 190801 (2016) . [32] Cosmo Lupo and Stefano Pirandola, “Ultimate Precision Bound of Quantum and Subwa velength Imaging, ” Physical Revie w Letters 117 , 190802 (2016) . [33] Mankei Tsang, “Conservati ve classical and quantum resolution limits for incoherent imaging, ” Journal of Modern Optics 65 , 1385–1391 (2018) . [34] Shan Zheng Ang, Ranjith Nair , and Mankei Tsang, “Quantum limit for two-dimensional resolution of two incoherent optical point sources, ” Physical Re view A 95 , 063847 (2017) . [35] Xiao-Ming Lu, Hari Krovi, Ranjith Nair, Saikat Guha, and Jeffre y H. Shapiro, “Quantum-optimal detection of one-v ersus- two incoherent optical sources with arbitrary separation, ” npj Quantum Information 4 , 64 (2018) . [36] J. ˇ Reh ´ a ˇ cek, M. Pa ´ ur , B. Stoklasa, Z. Hradil, and L. L. S ´ anchez-Soto, “Optimal measurements for resolution beyond the Rayleigh limit, ” Optics Letters 42 , 231–234 (2017) . [37] Fan Y ang, Ranjith Nair, Mankei Tsang, Christoph Simon, and Alexander I. Lvo vsky , “Fisher information for far-field linear optical superresolution via homodyne or heterodyne detection in a higher -order local oscillator mode, ” Physical Revie w A 96 , 063829 (2017) . [38] Ronan Kerviche, Saikat Guha, and Amit Ashok, “Fundamen- tal limit of resolving two point sources limited by an arbitrary point spread function, ” in 2017 IEEE International Symposium on Information Theory (ISIT) (IEEE, Aachen, Germany , 2017) pp. 441–445. [39] Andrzej Chrostowski, Rafał Demko wicz-Dobrza ´ nski, Marcin Jarzyna, and Konrad Banaszek, “On super-resolution imaging as a multiparameter estimation problem, ” International Journal of Quantum Information 15 , 1740005 (2017) . [40] J. ˇ Reh ´ a ˇ cek, Z. Hradil, B. Stoklasa, M. Pa ´ ur , J. Grov er , A. Krzic, and L. L. S ´ anchez-Soto, “Multiparameter quantum metrology of incoherent point sources: T ow ards realistic superresolution, ” Physical Re view A 96 , 062107 (2017) . [41] J. ˇ Reh ´ a ˇ cek, Z. Hradil, D. Koutn ´ y, J. Gro ver , A. Krzic, and L. L. S ´ anchez-Soto, “Optimal measurements for quantum spatial su- perresolution, ” Physical Re view A 98 , 012103 (2018) . [42] Mikael P . Backlund, Y oav Shechtman, and Ronald L. W alsworth, “Fundamental Precision Bounds for Three- Dimensional Optical Localization Microscopy with Poisson Statistics, ” Physical Re view Letters 121 , 023904 (2018) . [43] Carmine Napoli, Samanta Piano, Richard Leach, Gerardo Adesso, and T ommaso T ufarelli, “T ow ards Superresolution Surface Metrology: Quantum Estimation of Angular and Ax- ial Separations, ” Physical Re view Letters 122 , 140505 (2019) . [44] Zhixian Y u and Sudhakar Prasad, “Quantum Limited Super- resolution of an Incoherent Source Pair in Three Dimensions, ” Physical Re view Letters 121 , 180504 (2018) . [45] Sudhakar Prasad and Zhixian Y u, “Quantum-limited super- localization and superresolution of a source pair in three dimen- sions, ” Physical Re view A 99 , 022116 (2019) . [46] W alker Larson and Bahaa E. A. Saleh, “Resurgence of Rayleigh’ s curse in the presence of partial coherence, ” Optica 5 , 1382–1389 (2018) . [47] Mankei Tsang and Ranjith Nair , “Resurgence of Rayleigh’ s curse in the presence of partial coherence: comment, ” Optica 6 , 400–401 (2019) . [48] W alker Larson and Bahaa E. A. Saleh, “Resurgence of Rayleigh’ s curse in the presence of partial coherence: reply , ” Optica 6 , 402–403 (2019) . [49] Michael R. Grace and Saikat Guha, priv ate communication (2019). [50] Kent Bonsma-Fisher, W eng-Kian Tham, Hugo Ferretti, and Aephraim Steinberg, “Realistic sub-Rayleigh imaging with phase-sensitiv e measurements, ” arXiv:1906.00081 [physics, physics:quant-ph] (2019) . [51] Mankei Tsang, “Resolving starlight: a quantum perspective, ” arXiv:1906.02064 [astro-ph, physics:physics, physics:quant- ph] (2019) . [52] Zong Sheng T ang, Kadir Durak, and Alexander Ling, “Fault- tolerant and finite-error localization for point emitters within the diffraction limit, ” Optics Express 24 , 22004 (2016) . [53] W eng-Kian Tham, Hugo Ferretti, and Aephraim M. Steinberg, “Beating Rayleigh’ s Curse by Imaging Using Phase Informa- tion, ” Physical Re view Letters 118 , 070801 (2017) . [54] Martin Pa ´ ur , Bohumil Stoklasa, Zden ˇ ek Hradil, Luis L. S ´ anchez-Soto, and Jaroslav ˇ Reh ´ a ˇ cek, “ Achieving the ultimate optical resolution, ” Optica 3 , 1144 (2016) . [55] Fan Y ang, Arina T ashchilina, E. S. Moiseev , Christoph Simon, and A. I. Lv ovsk y , “Far-field linear optical superresolution via heterodyne detection in a higher-order local oscillator mode, ” Optica 3 , 1148 (2016) . [56] J. M. Donohue, V . Ansari, J. ˇ Reh ´ a ˇ cek, Z. Hradil, B. Stoklasa, M. Pa ´ ur , L. L. S ´ anchez-Soto, and C. Silberhorn, “Quantum- Limited T ime-Frequency Estimation through Mode-Selective Photon Measurement, ” Physical Re view Letters 121 , 090501 (2018) . [57] J. Hassett, T . Malhorta, M. A. Alonso, R. W . Boyd, S. M. Hashemi Rafsanjani, and A. N. V amiv akas, “Sub- Rayleigh Limit Localization with a Spatial Mode Analyzer, ” in F r ontiers in Optics / Laser Science (2018), paper JW4A.124 (Optical Society of America, 2018) p. JW4A.124. [58] Y iyu Zhou, Jing Y ang, Jeremy D. Hassett, Seyed Moham- mad Hashemi Rafsanjani, Mohammad Mirhosseini, A. Nick V ami vakas, Andrew N. Jordan, Zhimin Shi, and Robert W . Boyd, “Quantum-limited estimation of the axial separation of two incoherent point sources, ” Optica 6 , 534–541 (2019) . 15 [59] Erhan C ¸ ınlar , Pr obability and Stochastics (Springer , Ne w Y ork, 2011). [60] Donald L. Snyder and Michael I. Miller, Random P oint Pr o- cesses in T ime and Space (Springer-V erlag, New Y ork, 1991). [61] Michael Reed and Barry Simon, Methods of Modern Mathe- matical Physics. I: Functional Analysis (Academic Press, San Diego, 1980). [62] Mankei Tsang, “Quantum nonlocality in weak-thermal-light in- terferometry , ” Physical Re view Letters 107 , 270402 (2011) . [63] Alexander Meister , “Support estimation via moment estimation in presence of noise, ” Statistics 40 , 259–275 (2006) . [64] Charles F . Dunkl and Y uan Xu, Orthogonal P olynomials of Sev- eral V ariables (Cambridge Univ ersity Press, Cambridge, 2001). [65] D. J. Daley and D. V ere-Jones, An Intr oduction to the Theory of P oint Pr ocesses: V olume I: Elementary Theory and Methods (Springer-V erlag, New Y ork, 2003). [66] F . W . J. Olver , D. W . Lozier, R. F . Boisvert, and C. W . Clark, eds., NIST Handbook of Mathematical Functions (NIST and Cambridge Univ ersity Press, Cambridge, 2010). [67] Lokenath Debnath and Firdous Ahmad Shah, W avelet T rans- forms and Their Applications (Springer , New Y ork, 2015). [68] I. A. Ibragimov and R. Z. Has’minskii, Statistical Estimation: Asymptotic Theory (Springer , New Y ork, 1981). [69] Michał Parniak, Sebastian Bor ´ owka, Kajetan Boroszko, W oj- ciech W asilewski, Konrad Banaszek, and Rafał Demko wicz- Dobrza ´ nski, “Beating the Rayleigh Limit Using T wo-Photon Interference, ” Physical Re view Letters 121 , 250503 (2018) . [70] Martin Pa ´ ur , Bohumil Stoklasa, Jai Grover , Andrej Krzic, Luis L. S ´ anchez-Soto, Zden ˇ ek Hradil, and Jarosla v ˇ Reh ´ a ˇ cek, “T empering Rayleigh’ s curse with PSF shaping, ” Optica 5 , 1177–1180 (2018) . [71] M. Pa ´ ur , B. Stoklasa, D. Koutn ´ y, J. ˇ Reh ´ a ˇ cek, Z. Hradil, J. Grover , A. Krzic, and L. L. S ´ anchez-Soto, “Reading out Fisher information from the zeros of the point spread function, ” Optics Letters 44 , 3114–3117 (2019) . [72] Carl W . Helstrom, Quantum Detection and Estimation Theory (Academic Press, New Y ork, 1976). [73] Alexander S. Holev o, Pr obabilistic and Statistical Aspects of Quantum Theory (Edizioni della Normale, Pisa, Italy , 2011). [74] Akio Fujiwara, “Geometry of Quantum Estimation Theory , ” in Asymptotic Theory of Quantum Statistical Infer ence: Selected P apers , edited by Masahito Hayashi (W orld Scientific, Singa- pore, 2005) Chap. 18, pp. 229–286. [75] Konrad Schm ¨ udgen, The Moment Pr oblem (Springer, Cham, 2017). [76] J. D. Gorman and A. O. Hero, “Lower bounds for parametric estimation with constraints, ” IEEE T ransactions on Information Theory 36 , 1285–1301 (1990) . [77] Roger A. Horn and Charles R. Johnson, Matrix Analysis (Cam- bridge Univ ersity Press, Cambridge, 1985). [78] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Riv est, and Clifford Stein, Introduction to Algorithms (MIT Press, Cambridge, 2009).