MANAGING ADAPTIVE MEASUREMENT FOR HIGH-RESOLUTION MEASUREMENT

Abstract

Imaging a distribution of one or more optical sources includes: receiving respective optical signals from a spatial mode sorter during each of two or more detection intervals of time; after each of the two or more detection intervals of time, processing information based at least in part on: (1) the respective optical signal received in the corresponding detection interval of time, and (2) a first set of models comprising a set of distributions related to one or more optical sources, each model corresponding to a different number of optical sources in the distribution, and configuring the spatial mode sorter based at least in part on the processing; and providing an estimated measurement characterizing the distribution of one or more optical sources based at least in part on the processed information. The processing after at least one of the two or more detection intervals of time includes computing an eigen-projection.

Description

TECHNICAL FIELD

This disclosure relates to managing adaptive measurement for high-resolution measurement.

BACKGROUND

Some passive imaging architectures employ a digital focal plane followed by electronic-domain post-processing. However, this approach can be inefficient, especially for resolving features with extent smaller than the Rayleigh diffraction limit. Using quantum information theory, suitable spatial modal measurement(s) can significantly outperform the traditional diffraction-limited imaging for estimating two sub-Rayleigh spaced incoherent point sources (e.g., below the diffraction limit).

SUMMARY

In one aspect, in general, a method for imaging a distribution of one or more optical sources includes: receiving respective optical signals from a spatial mode sorter during each of two or more detection intervals of time; after each of the two or more detection intervals of time, processing information based at least in part on: (1) the respective optical signal received in the corresponding detection interval of time, and (2) a first set of models comprising a set of distributions related to one or more optical sources, each model corresponding to a different number of optical sources in the distribution, and configuring the spatial mode sorter based at least in part on the processing; and providing an estimated measurement characterizing the distribution of one or more optical sources based at least in part on the processed information. The processing after at least one of the two or more detection intervals of time includes computing an eigen-projection.

Aspects can include one or more of the following features.

The configuring after each of the two or more detection intervals of time configures the spatial mode sorter to project the respective optical signals onto a basis of the computed eigen-projection.

The first set of models includes a set of spatial distributions for the optical sources.

The first set of models further includes a set of Bayesian prior probability distributions for the set of spatial distributions.

The set of Bayesian prior probability distributions for the set of spatial distributions includes a set of Gaussian distributions.

A set of hyper-parameters for the Gaussian distributions are based at least in part on a result of an expectation maximization calculation that is based at least in part on the processed information.

The first set of models includes a set of brightness distributions for the optical sources.

The first set of models further includes a set of Bayesian prior probability distributions for the set of brightness distributions.

The set of Bayesian prior probability distributions for the set of brightness distributions includes a set of Dirichlet distributions.

A set of hyper-parameters for the Dirichlet distribution is based at least in part on a result of an expectation maximization calculation that is based at least in part on the processed information.

The eigen-projection is the eigenvectors of a symmetric logarithmic derivative operator.

The symmetric logarithmic derivative operator is based at least in part on a set of one or more operators constructed from a single-parameter inference setting.

The eigen-projection is the eigenvectors of an operator constructed from a Bayesian inference setting.

The operator constructed front a Bayesian inference setting is based at least in part on a set of one or more operators constructed from a single-parameter Bayesian inference setting.

The eigen-projection comprises a Personick eigen-projection.

The processed information includes a second set of models.

The second set of models comprising a second set of distributions related to one or more optical sources, determined by (1) the first set of distributions related to one or more optical sources and (2) respective optical signals received in a previous detection interval of time.

The processing after at least one of the two or more detection intervals of time includes computing a quantum Fisher information matrix associated with the respective optical signals.

The processing after at least one of the two or more detection intervals of time includes computing a modified quantum Fisher information matrix derived in a Bayesian inference setting and associated with the respective optical signals.

In another aspect, in general, one or more non-transitory computer-readable media, having instructions stored thereon that, when executed by a computer system, cause the computer system to perform operations including: receiving respective optical signals from a spatial mode sorter during each of two or more detection intervals of time; after each of the two or more detection intervals of time, processing information based at least in part on: (1) the respective optical signal received in the corresponding detection interval of time, and (2) a first set of models comprising set of distributions related to one or more optical sources, each model corresponding to a different number of optical sources in the distribution, and configuring the spatial mode sorter based at least in part on the processing; and providing an estimated measurement characterizing the distribution of one or more optical sources based at least in part on the processed information. The processing after at least one of the two or more detection intervals of time includes computing an eigen-projection.

In another aspect, in general, an apparatus for imaging a distribution of one or more optical sources includes: a spatial mode sorter that defines a configurable basis comprising a set of spatial modes onto which an coming optical signal is projected; and a control module. The control module is configured to: receive respective optical signals from the spatial mode sorter during each of two or more detection intervals of time; after each of the two or more detection intervals of time, process information based at least in part on: (1) the respective optical signal received in the corresponding detection interval of time, and (2) a first set of models comprising a set of distributions related to one or more optical sources, each mod& corresponding to a different number of optical sources in the distribution, and configure the spatial mode sorter based at least in part on the processing; and provide an estimated measurement characterizing the distribution of one or more optical sources based at least in part on the processed information. The processing after at least one of the two or more detection intervals of time includes computing an eigen-projection.

Aspects can have one or more of the following advantages.

Some of the techniques described herein use an eigen-projection adaptive algorithm for super-resolution imaging. For example, an adaptive measurement scheme can resolve point-source constellations with sub-Rayleigh separations by employing certain modal measurements. The measurements can be performed using classical measurement techniques, but may be based on certain aspects of quantum estimation theory (e.g., quantum-inspired optical super-resolution adaptive imaging). These measurement techniques can achieve increased performance, as described in more detail below.

Other features and advantages will become apparent from the following description, and from he figures and claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The disclosure is best understood from the following detailed description when read in conjunction with the accompanying drawing. It is emphasized that, according to common practice, the various features of the drawing are not to-scale. On the contrary, the dimensions of the various features are arbitrarily expanded or educed for clarity.

FIG. 1 is a schematic diagram of an example adaptive measurement system.

FIG. 2 is a schematic diagram of an example adaptive measurement procedure.

FIG. 3A and FIG. 3B are plots of example point source locations and a pseudo-imaging metric, respectively.

FIG. 4 is a schematic diagram of an example SLD approach.

FIG. 5A is a plot of an example of estimation given by direct imaging and Personick projection measurement.

FIG. 5B is a plot of the distribution of the point emitter position errors obtained direct imaging versus Personick projection, when the number of emitters is known exactly.

FIG. 5C is a plot of the distribution of the number of point emitters estimated by direct imaging and Personick projection.

FIG. 5D is a plot of the distribution of the point emitter position errors obtained with direct imaging versus Personick projection, when the number of emitters is not previously known.

FIG. 6A is a schematic diagram of an implementation of a spatial mode sorter system.

FIG. 6B is a schematic diagram illustrating the sorting of a first spatial mode.

FIG. 6C is a schematic diagram illustrating the sorting of a second spatial mode.

FIG. 6D is a schematic diagram illustrating the sorting of a third spatial mode.

FIG. 7 is a schematic diagram of a second implementation of a spatial mode sorter.

FIG. 8 is a schematic diagram of a third implementation of a spatial mode sorter.

DETAILED DESCRIPTION

A measurement method is designed to enhance the performance of estimating the brightness and location of incoherent point sources, including resolving the point pairs which are in the sub-Rayleigh regime, in a constellation through the light collected by incoherent imaging systems like microscopes or telescopes. In some implementations, the design is based on the fact that for a single parameter estimation problem, the eigen-projection measurement of the symmetric logarithmic derivative (SLD) operator achieves the quantum Fisher information (QFI) limit. For the 3N parameters (N for brightness and 2N for the location on the 2D focal plane) of a N-point constellation, the eigenvector with the largest eigenvalue of the current state of estimation is identified as the single parameter and a set SLD eigen-projection is constructed for the measurement of the next stage. Other implementations are based on another type of eigen-projection measurement, as described in more detail below.

One example of the algorithm includes two main stages: initialization, and SLD eigen-projection. In this example, in the initialization stage, direct imaging measurement with an expectation maximization (EM) algorithm is used to get a set of pre-estimation for the 3N parameters. Based on the pre-estimation, a set of density operators can be created using the grouping algorithm and the density operators are carried to the next stage of measurement. The initialization stage can be performed as described herein, or by any technique that provides a suitable pre-estimation. In the SLD eigen-projection stage, the adaptive measurement is divided into different cycles. In each cycle, a density operator is selected from the pool of models and the QFI matrix of that operator is calculated. The eigenvector with the largest eigenvalue is selected as the single parameter which is used for the construction of the SLD eigen-projection measurement. A certain amount of photons from the incoherent sources are measured in the spatial modes given by the eigen-projection of SLD. The density operators are updated under the Bayesian framework and one of them is selected according to their likelihood for the next cycle. When photons are exhausted, the parameters in the density operator being picked is used as the final estimation. Other techniques can be used to select the density operator and the parameter. The selection schemes may be task specific.

Referring to FIG. 1, an example of a measurement system 100 includes an adaptive receiving system 102 that includes an adaptive engine 104 (e.g., including an optical front-end and a processing module implemented on a special-purpose or general-purpose processor) for receiving an optical input 103 and producing measurement information 105 as output. The adaptive engine 104 receives model information 106 defining multiple models for a potential distribution of optical sources to be imaged. For example, a set of models 108 are stored in a storage system 110. Between a series of detection intervals of time, the adaptive engine 104 configures a configurable spatial mode sorter 112 to provide optical signals that change based at least in part on the adaptive computations being performed, as described in more detail herein. Examples of some aspects of this part of the procedure (e.g., spatial-mode sorting) are described in more detail below.

Without intending to be bound by theory, simulation has been done to compare the performance of this scheme against a traditional diffraction limited direct imaging under different metrics. The result shows that for the Jaccard index metrics, which calculate the Jaccard index for the point sources together with RMSE for the parameters, and the pseudo-imaging metrics, which calculate the overlapping integral of the image given by the ground truth of the scene and the estimated scene through a pseudo-aperture, the SLD eigen-projection outperforms the traditional method.

FIG. 2 shows an example of an adaptive measurement procedure 200. Given a constellation composed of point-sources as shown in the upper left of FIG. 2, an example of imaging problem can be stated as estimating: 1) the number of point sources, 2) their relative brightness and 3) positions, given prior information on the maximum number of point sources P_max. In some implementations, a goal is to resolve all point sources at all separations including those in sub-Rayleigh regime (super-resolution problem). We employ a diffraction-limited lens with a Gaussian aperture, and propose an adaptive modal projection measurement scheme. The adaptive scheme includes three stages: initial measurement (initialization), adaptive measurements and termination.

At a high level, the adaptive measurement procedure 200 images a constellation of an unknown number of point sources acid constructs a model for the location of the point sources, with one model for each possible number of point sources in the constellation. The number of point sources may be bounded by prior knowledge. A projective measurement is constructed from the model and the constellation is further imaged. The additional information from the imaging is then used to update each model. This process can be iterated multiple times until a goal is reached (e.g., until no further imaging is allowed), after which the model with the highest likelihood is used and an estimate for the number of point sources and their locations is determined.

In the initialization stage, we allocate a pre-defined number of photons to a traditional imaging measurement. An initial estimate of the constellation parameters is obtained using this imaging measurement using the expectation maximization (EM) algorithm. Based on the estimated parameters and a grouping algorithm we have developed, we construct a set of models of the constellation with varying number of point sources ranging from NEM, number of point sources estimated by the EM algorithm, to P_max. Each model is represented by a density operator with incoherent mixture of spatial coherent states. For each model, a prior distribution is assigned to each position and brightness parameter.

Modal projection measurements are performed in the next stage. In some implementations, the design of the modal basis is based on the fact that for a single parameter problem, the eigen-projection of the symmetric logarithmic derivative (SLD) operator achieves the quantum limit. However, the problem here is a multi-parameter problem and furthermore the SLD eigen-projection is parameter dependent. Thus, the system first selects a model with the highest likelihood given the photons collected so far. Next, the system calculates the quantum Fisher information (QFI) matrix for that model (parameters) and pick the eigenvector with highest eigenvalue as the hypothetical parameter γ that is used to construct the SLD operator. Measuring the photons using the SLD eigen-projection of γ, the posterior is updated and serves as the prior for the next measurement. The measurement scheme terminated when the photon/exposure time budget is exhausted and the model with highest likelihood is selected.

Referring to FIG. 3A, we use an example with 8 point sources to demonstrate the performance of an example adaptive measurement scheme. The point sources are generated (uniformly) within a 3 RL (Rayleigh length)×3 RL square field of view with relative brightnesses following a uniform distribution. The maximum number of point sources P_max is set 10 and the total photon budget is 10000. The SLD eigen-projection measurement is initialized by using 500 photons on an imaging measurement followed by EM algorithm to estimate the initial parameters. The remaining 9500 photons are detected by the sequentially adaptive SLD eigen-projection measurements. For some imaging implementations, the system uses an EM algorithm for estimating the constellation. We observe that for sub-Rayleigh separated point sources (dashed box), the SLD eigen-projection measurements are able to resolve the sources while the imaging aggregates them as a single source.

FIG. 3B shows the resulting imaging performance evaluation using a metric/score M_x defined as:

$M_{\times }\left(\sigma \right)=\int \sqrt{I\left(\theta ;\sigma \right)I\left(\hat{\theta };\sigma \right)}d\overrightarrow{r},$

where θ is the true parameter value, {circumflex over (θ)} is the estimated parameter value and σ is the width of the pseudo-PSF.

This imaging metric measures the overlapping integral of the pseudo-intensity distribution I(θ, σ) for true and estimated parameters. In this example, the intensity distribution can be defined as super-position of intensity of point sources blurred by a Gaussian PSF (with width σ). As σ decreases the imaging metric becomes more sensitive, while its overall magnitude also decreases. In this prophetic example, the system averages around 900 constellations at different σ, where P_max=10 and the true number of point sources P ranging from 5 to 10. We can see that the SLD eigen-projection measurement (upper line) outperforms the imaging (lower line) for all σ.

These results show that the proposed adaptive imaging scheme is capable of optical super-resolution even in the noise-limited regime. Implementations of such a scheme have potential applications in space surveillance, astronomy, and biological imaging (fluorescent microscopy).

A variety of implementations of the approach described herein can be used. FIG. 4 shows a high level summary 400 of some such implementations. Other implementations are also possible.

In Bayesian estimation theory the estimator {circumflex over (θ)}=E[θ|l] attains the minimum mean squared error (MMSE) for estimating a scalar parameter of interest θ from the observation of l through a noisy channel P_l|θ, given a prior P_θ on θ. In quantum, e.g., optical and atomic, sensing tasks the user gets ρ_θ, the quantum state that encodes θ. The user chooses a measurement, a positive-operator valued measure (POVM)Π_l, which induces the channel P_l|θ=Tr(ρ_θΠ_l) to the measurement outcome l, on which the aforesaid classical MMSE estimator is employed. There exists an optimum POVM Π_l that minimizes the MMSE over all possible measurements, and the corresponding MMSE. This result is an alternative to the quantum Fisher information (QFI), which lower bounds the variance of an unbiased estimator over all measurements, when P_θ is unavailable. For multi-parameter estimation, i.e., when θ is a vector, in Fisher quantum estimation theory, the inverse of the QFI matrix provides an operator lower bound to the covariance of an unbiased estimator. However, there has been little work on quantifying quantum limits and measurement designs for multi-parameter quantum estimation in the Bayesian setting.

Some implementations build upon this result to construct a Bayesian adaptive measurement scheme for multi-parameter estimation when P copies of ρ_θ are available. We illustrate an application to localizing a cluster of point emitters in a highly sub-Rayleigh angular field-of-view, an important problem in fluorescence microscopy and astronomy. An example algorithm translates to a multi-spatial-mode transformation prior to a photon-detection array, with electro-optic feedback to adapt the mode sorter. We show that this receiver outperforms quantum-noise-limited focal-plane direct imaging.

Without intending to be bound by theory, an example of a sensing and imaging measurement design problem can be abstracted to having access to an P-copy quantum state ρ(θ)^⊗N encoding M parameters of interest, denoted by the vector θ=[θ₁, θ₂, . . . , θ_M]^T. The receiver realizes a quantum measurement, described by a positive-operator-valued measure (POVM) {Π_l}, operating on a single copy of ρ(θ), resulting in a vector-valued measurement outcome l=[l₁, l₂, . . . , l_N]^T. The choice of a POVM induces the (classical') measurement channel defined by the probability density p(l|θ)=Tr(ρ(θ)Π_l). When more than one copy of the quantum state is available (N≥2), the receiver can: (1) in the most general setting choose a joint-measurement POVM {Π_l(N)} acting collectively on ρ(θ)^⊗N, producing the outcome l_(N): (2) employ the so-called local operations and classical communications (LOCC) scheme, where each copy of the state is measured by an independent measurement, where the POVM {Π_l^τ} acting on the τ^th copy of ρ(θ) is chosen based on the information available from the previous set of measurement outcomes {l⁽¹⁾, l⁽²⁾, . . . , l^(τ−1)}, 1≤τ≤N; or (3) use independent identical measurements on each copy of the state, described by the POVM {Π_l}. No matter which strategy the receiver may use, after measuring all P copies, the receiver generates an estimate of θ, i.e., {circumflex over (θ)} (l_set) where l_set=l_(N) for case (1) above, and l_set=[l⁽¹⁾, l⁽²⁾, . . . , l^(N)] for cases (2) and (3) above. The receiver chooses the estimator to optimize a desired objective/loss function. A natural choice of the objective function associated with sensing and imaging tasks is mean (expected) squared-error (MSE), E[∥θ−{circumflex over (θ)}∥²].

For any given measurement POVM {Π_l}, assuming strategy (3) above, i.e., the same measurement acts on each copy of ρ(θ), the problem reduces to the standard classical estimation theory problem of estimating θ from P i.i.d. samples of l, each produced by p(l|θ). The covariance Cov({circumflex over (θ)}(l_set), θ) for any unbiased estimator {circumflex over (θ)} of θ is lower bounded by Σ_C. This means Cov({circumflex over (θ)}(l_set), θ)−Σ_C is a semi-positive definite matrix, denoted compactly as Cov({circumflex over (θ)}(l_set),θ)≥Σ_C. The receiver's task is to pick the optimal estimator {circumflex over (θ)}^opt (l_set) on the measurement l_set, such that Cov({circumflex over (θ)}^opt (l_set), θ) saturates the bound Σ_C when permissible.

Tools of quantum estimation theory let us find a tighter lower bound to Cov({circumflex over (θ)}(l), θ), which automatically optimizes over all physically-permissible choices of a POVM {Π_l} (again, assuming that the same measurement is used to detect each copy of ρ(θ)). The Cov ({circumflex over (θ)}(l), θ) is lower bounded by Σ_Q (a quantum bound.), which itself is an infimum of all bounds Σ_C associated with all possible choices of {Π_l}. For certain cases (for example when θ is a single scalar parameter), quantum estimation theory also tells us the optimal receiver POVM {Π_l^(opt)}. Once the optimal receiver is chosen, it uses the optimal estimator {circumflex over (θ)}^opt (l_set) using standard classical estimation tools, such that Cov ({circumflex over (θ)}^opt (l_set), θ) saturates Σ_Q when permissible. In general Cov({circumflex over (θ)}(l_set), θ)≥Σ_C≥Σ_Q, where Σ_C could correspond to any POVM choice.

The aforementioned lower bounds on the covariance of multi-parameter estimators can be defined within the statistical inference frameworks the frequentist, i.e., Fisherian (with no prior), or the Bayesian (with prior p(θ)) inference settings. We review below some known bounds in both the settings.

In the Fisherian (frequentist) framework, when no prior p(θ) is available, the Cramer-Rao lower bound (CRLB) Σ_C on the covariance Cov({circumflex over (θ)} (l), θ) of an unbiased estimator is given by the inverse of the Fisher information (FI) matrix I.

$\begin{array}{cc}{I}_{\mathrm{ij}}=\int \left[\frac{\partial }{\partial {\theta }_i}\ln p\left(l❘\theta \right)\right]\left[\frac{\partial }{\partial {\theta }_j}\ln p\left(l❘\theta \right)\right]p\left(l❘\theta \right)\mathrm{dl},& \left(1\right)\end{array}$

with I≤(i,j)≤M, and the likelihood p(l/θ)=Tr(ρ(θ)Π_l). The quantum version of this lower bound Σ_Q, which only depends on ρ(θ) (since the measurement Π_l is automatically optimized over) is given by the inverse of the quantum Fisher information (QFI) matrix Q, the elements ofwhich are:

$\begin{array}{cc}{Q}_{\mathrm{ij}}=\frac{1}{2}\mathrm{Tr}\left[\rho \left(\theta \right)\frac{L_i{L}_j+L_j{L}_i}{2}\right],& \left(2\right)\end{array}$

where L_i is the symmetric logarithmic derivative (SLD) operator which can be evaluated from the implicit relationship:

$\begin{array}{cc}2\frac{\partial }{\partial {\theta }_i}\rho \left(\theta \right)=\rho \left(\theta \right)L_i+L_i\rho \left(\theta \right),& \left(3\right)\end{array}$

with 1≤i≤M. Thus, we have Cov({circumflex over (θ)} (l_set), θ)≥I⁻¹≥Q⁻¹ in the Fisher framework. For N-copy i.i.d. measurement of ρ(θ)^⊗N, both the classical and quantum bounds above are lowered by a factor of 1/N. The classical one is asymptotically attained by the maximum likelihood estimator (NILE). The quantum CRLB (Q⁻¹) is not saturable in general for M>1.

Given a prior p(θ) on the pan Teter vector θ, the Bayesian Cramer-Rao lower bound (BCRLB) Σ_C is given by:

$\begin{array}{cc}{\sum }_C=\int p\left(\theta \right){\mathrm{\theta \theta }}^{T}d\theta -J,& \left(4\right)\end{array}$

where the M-by-M matrix J is defined as:

$\begin{array}{cc}{J}_{\mathrm{ij}}=\int \frac{\left[\int p\left(l,\theta \right){\theta }_{i}d\theta \right]\left[\int p\left(l,\theta \right){\theta }_{j}d\theta \right]}{p\left(l\right)}\mathrm{dl},& \left(5\right)\end{array}$

and p(l, θ)=p(l|θ)p(θ) is the joint distribution of l and θ. For the quantum version of this lower bound, we first define the following operators, for 1≤i≤M and k=0,1,2:

$\begin{array}{cc}{\Gamma }_{i,k}=\int d\theta p\left(\theta \right)\rho \left(\theta \right){\theta }_i^k,& \left(6\right)\end{array}$

and operators B_i, l≤i≤M, that satisfy:

$\begin{array}{cc}2{\Gamma }_{i,1}={\Gamma }_0{B}_i+B_i{\Gamma }_0.& \left(7\right)\end{array}$

For k=0, Γ_i,0=Γ_j,0, ∀(i,j), thus we can drop the first index and denote it as Γ₀=∫ dθp(θ)ρ(θ), the average received state. The quantum BCRLB Σ_Q can be written as:

$\begin{array}{cc}{\sum }_Q=\int p\left(\theta \right){\mathrm{\theta \theta }}^{T}d\theta -G,& \left(8\right)\end{array}$

where

$\begin{array}{cc}{G}_{\mathrm{ij}}=\frac{1}{2}\mathrm{Tr}\left[{\Gamma }_0\frac{B_i{B}_j+B_j{B}_i}{2}\right].& \left(9\right)\end{array}$

Thus in a Bayesian inference framework, we have Cov({circumflex over (θ)} (l),θ)≥Σ_C≥Σ_Q.

To achieve the quantum bound, an optimal measurement is required (i.e. an optimal choice POVM, that acts on each copy of ρ(θ)). For a single parameter problem (M−1), the projective measurement onto the eigenvectors of the SLD operator L in Eq. (3) saturates the Fisher bound, i.e., the Σ_C for the SLD measurement equals Σ_Q. Likewise, the Bayesian quantum bound on the covariance is saturated, for the case of a single parameter (M−1) by a projective measurement onto the eigenvectors of the operator B in Eq. (7).

For multi-parameter estimation, if the operators associated with parameter θ_i: L_i and B_i, 1≤i≤M commute with one another, for the Fisher and Bayesian frameworks respectively, the corresponding covariance bound can be saturated by the above-said measurements, calculated by evaluating the eigenvectors of L_i or B_i, respectively (which i does not matter as they are simultaneously diagonal). However, if the operators do not commute, which is the case in general, a measurement that is jointly optimal for all parameters may not exist and is likely to be challenging to derive.

In the quantum case, the Holevo Cramer-Rao bound (HCRB) is the most fundamental scalar lower bound on the weighted mean square error Tr(WCov({circumflex over (θ)}, θ)), for a positive definite W. The HCRB represents the best precision attainable with a collective measurement (discussed as case (1) above) on an asymptotically large number of identical copies of ρ(θ).

We propose a sequential adaptive measurement scheme (an LOCC) within a full Bayesian inference framework. We leverage tools from Bayesian quantum estimation theory. In one example, we employ our measurement scheme to the problem of localizing an unknown number of point-emitters placed in a sub-Rayleigh (below diffraction-limit) field of view in an optical imaging context. This application is motivated by the fact that traditional direct focal-plane imaging, which employs intensity measurements followed by electronic-domain processing, is known to be highly sub-optimal in the sub-Rayleigh regime. We compare our quantum-inspired adaptive sequential measurement design with the direct imaging technique to quantify the significant performance improvement in optical resolution achieved by our scheme,

Consider a system in the state described by the density operator:

$\begin{array}{cc}{\rho }^{\prime }\left(\theta \right)=\sum _{i=1}^P{b}_i\left(\theta \right)❘{\psi }_i\left(\theta \right)\rangle \langle {\psi }_i\left(\theta \right)❘,& \left(10\right)\end{array}$

where θ=[θ₁, θ₂, . . . , θ_M are the parameters of interest, |Ψ_i(θ) and b_i(θ) are the parameter-dependent pure states and the corresponding weights respectively. P itself, in general, is an unknown parameter (positive integer) such that: P_min≤P≤P_max. Here we assume that P is upper bounded by P_max, i.e., a prior on P. If the lower bound P_min is not known/available, we can set it to 1. When P_min≠P_max, both P and θ parameters need to be estimated. On the contrary, if P_min=P=P_max, i.e., P is known a priori exactly and we only need estimate the parameters θ.

A single copy of a quantum state ρ′(θ), defined in Eq. (10), can be measured with a POVM {Π_l} such that the probability of getting a scalar outcome l is p(l|θ)=Tr[ρ′(θ)Π_l]. Now, let K denote the block-length (of copies of ρ′(θ)) over which the measurement Π_l stays the same. The density operator ρ(θ)=ρ′(θ)^⊗K and the probability of observing the measurement outcome l=[l₁, l₂, . . . , l_K]^T is p(l|θ) =Tr[ρ(θ)Π_l]=Π_i=1^K Tr[ρ′(θ) Π_li], where Π_l Π_l1 ⊗ . . . ⊗Π_lK.

To illustrate the proposed scheme, we begin with the P known exactly case. An extension where we relax this prior on P is discussed in the next section. Let us take N=K×S to be the total number of copies of ρ′(θ) available to us. We adapt the measurement choice S times, denoted by τ as the iteration index, 1≤τ≤S. Thus, at the τ^th iteration of the sequential measurement, a single copy of ρ(θ) is measured with the POVM {Π_l^(τ} yielding the outcome vector l^(τ)=[l₁^(τ), l₂^(τ), . . . , l_K^(τ)]^T. At the end of the sequential measurement scheme, a S-copy state ρ(θ)^⊗S has been measured. The parameter estimate {circumflex over (θ)} ^(τ), after the τ^th sequential measurement is denoted by {circumflex over (θ)} ^(τ)=[{circumflex over (θ)} ₁^(τ), {circumflex over (θ)} ₂^(τ), . . . , {circumflex over (θ)} _M^(τ)]^T. In a Bayesian inference setting, the parameter estimate {circumflex over (θ)} ^(τ) 0 is given by posterior mean: {circumflex over (θ)} ^(τ)=E_θ(τ) [p(θ^(τ)|l^(τ))]. The posterior p(θ^(τ)|l^(τ))=p(l^(τ)|θ^τ−1))·p(θ^(τ−1))/p(l^(τ)), where p(θ|l^(τ−1)) at the previous τ^th iteration. Note that this prior p(θ^τ−1)) in turn equals the posterior p(θ)|l^(τ−1)) at the previous τ−l^th iteration, terminating at p(θ⁽¹⁾)=p(θ). The density operator at the τ^th iteration is represented as ρ(θ^(τ)). Now what remains to be determined is how we choose the POVM {Π_l^(τ)} at the τ^th iteration.

It is known that for a single parameter estimation problem, the eigen-projection measurement of B₁ in Eq. (7) saturates the quantum bound Σ_Q, which reduces to a lower bound of the variance of the scalar parameter. The minimum mean square error (MMSE) is given by Σ_Q=Tr[Γ_1,2−B₁Γ_1,1], where Γ_i,j are defined in Eq. (6). We refer to this measurement as the Personick projection in this work. For the multi-parameter problem, the counter-part of Tr[B₁Γ_1,1] is a matrix G in Eq. (9). If all B_i operators commute, the quantum optimal measurement is given by the eigen-projections onto any of the B_i operators, however, there is no such guarantee that the optimal measurement for all parameters exists or can be found. At the τ^th iteration of sequential measurement we propose to use a single parameter γ^(τ), defined as a linear combination of the M parameters given by the the eigenvector of the matrix G with the smallest eigenvalue. Note that this matrix G is defined per Eq. (9) for the density operator ρ({circumflex over (θ)} ^(τ−1)). The scalar parameter γ^(τ) is used to construct the operator B_γ(τ). The corresponding Personick projection constructed using B_γ(τ) is chosen as the POVM {Π_l^(τ} at the τ^th iteration. The sequential measurements are terminated when we have used all the P available copies of ρ′(θ).

If the scalar P in Eq. (10) is unknown, we can employ and initialize multiple models of density operators ρ(θ_P) with the corresponding prior p(θ_P;), where θ_P=[θ₁, θ₂, . . . , θ_MP]^T for P_min≤P≤P_max. In such a scenario, the number of parameters M_P for each model corresponding to a P can be different in general. In τ^th iteration of the sequential measurement, one mod& is selected and used to construct the Personick measurement. After model selection, the measurement scheme defined in the previous section can be applied. Note that at τ^th iteration, not only selected model but all the models are updated using the measurement outcome l^(τ). When the sequential measurements eventually terminate, we pick a model using same model selection criteria described above and compute the final multi-parameter estimate as the posterior mean.

We apply an adaptive sequential measurement scheme to estimate the location d relative brightness of incoherent point emitters in a cluster/constellation. The quantum state of photons incident on the image plane through a simple lens is given by the density operator ρ_full:

$\begin{array}{cc}{\rho }_{\mathrm{full}}=\left(1-\epsilon \right)❘0\rangle \langle 0❘+{\mathrm{\epsilon \rho }}^{\prime }+O\left({\epsilon }^2\right),& \left(11\right)\end{array}$

where |0 is the vacuum state, ρ′ is the single photon state density operator, which has the for .1 of Eq. (10), and ϵ is the total number of photons arriving on the image plane within the coherence time of the source. Assuming that ϵ≤≤1 (valid for weak thermal source), the photon states with order O(ϵ²) are negligible. As the vacuum state |0 provides no information, we can focus on ρ′. Thus, the components of Eq. (10) have the following meaning: P is the number of emitters, {b_i}_i=1^P are the relative brightness of each source (sum normalized to 1) and the states |Ω_i are given by:

$\begin{array}{cc}❘{\psi }_i\rangle ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\psi \left(x-x_i,y-y_i\right)❘x,y\rangle \mathrm{dxdy},& \left(12\right)\end{array}$

such that (x_i, y_i) are the scaled coordinates of the i^th point source on e image plane. The point spread function (PSF) Ω(x, y) of the system is set to a 2D Gaussian function:

$\begin{array}{cc}\psi \left(x,y\right)=\frac{1}{\sqrt{2{\mathrm{\pi \sigma }}_x{\sigma }_y}}\exp \left(-\frac{x^2}{4{\sigma }_x^2}-\frac{y^2}{4{\sigma }_y^2}\right),& \left(13\right)\end{array}$

where σ_x and σ_y are the standard deviation (a measure of width) of the PSF in x and y direction respectively. For a given PST, σ_x and σ_y are known parameters and set to σ_y=σ_y in our study. We define the full width at half maximum (proportional to σ_x) of the PSF as Rayleigh length (rl) in our analysis.

The parameters of interest in this problem are the position and relative brightness of the P emitters, i.e. θ=[x₁, . . . , x_P, y₁, . . . , y_P, b₁, . . . , b_P]^T=[x,y,b]^T. For the positions [x,y]^T, we use independent Gaussian distribution as prior:

$\begin{array}{cc}p\left(x,y\right)=\prod _i^P𝒩\left(x_i;{\stackrel{\_}{x}}_i,{\stackrel{\_}{\sigma }}_{x_i}\right)𝒩\left(y_i;{\stackrel{\_}{y}}_i,{\stackrel{\_}{\sigma }}_{y_i}\right),& \left(14\right)\end{array}$

where for 1≤i≤P, x_i, y_i, σ_xi, σ_yi are the mean and standard deviation of the position parameters x_i and y_i respectively. For the brightness b^T parameters a Dirichlet distribution is used as prior: p(b)=Dir(b; a), where a=[a₁, . . . , a_P]^T are the hyper-parameters of the Dirichlet distribution. Thus, the overall prior is expressed as: p(x,y,b)=p(x,y)p(b).

Now we have defined all relevant detail (i.e., photon state density operator, prior distribution) for the proposed adaptive sequential measurement scheme described in the previous section. As p(x , b) is not a conjugate prior for the Poisson likelihood, we update the hyper-parameters of the prior distribution at τ^th iteration to derive the posterior, which assumes the role of the prior in the next τ+l^th iteration. The prior hyper-parameters are: h=[x₁, . . . , x_P, y₁, . . . , y_P, σ_x1, . . . , σ_XP, σ_y1, . . . , σ_yP, a₁, . . . a_P, δ]^T=[x,y,σ_x, σ_y,a, δ]^T. Here, δ is another hyper-parameter associated with the brightness prior distribution explained later. To update the hyper-parameters of the position prior at the τ+l^th iteration, we use the first- and the second-moments of the posterior distribution at the τ^th iteration:

$\begin{array}{cc}{\stackrel{\_}{\alpha }}_i^{\left(\tau +1\right)}=\int {\alpha }_{i}p\left({\theta }^{\left(\tau \right)}❘l^{\left(\tau \right)};h^{\left(\tau \right)}\right)d{\theta }^{\left(\tau \right)},& \left(15\right)\end{array}$ $\begin{array}{cc}{\stackrel{\_}{\sigma }}_{{\alpha }_i^{\left(\tau +1\right)}}^2=\int {\left[{\alpha }_i-{\alpha }_i^{\left(\tau +1\right)}\right]}^{2}p\left({\theta }^{\left(\tau \right)}❘l^{\left(\tau \right)};h^{\left(\tau \right)}\right)d{\theta }^{\left(\tau \right)},& \left(16\right)\end{array}$

where α represents x or y coordinate.

For the hyper-parameters a^T of the brightness prior, expectation maximization approach is used. We first find the mean of the brightness vector as:

$\begin{array}{cc}{\hat{b}}_i^{\left(\tau +1\right)}=\int b_{i}p\left({\theta }^{\left(\tau \right)}❘l^{\left(\tau \right)};h^{\left(\tau \right)}\right)d{\theta }^{\left(\tau \right)}.& \left(17\right)\end{array}$

Then, α^T is updated such that {circumflex over (b)}^(τ+1) becomes the mode

$\begin{array}{cc}\begin{array}{c}{a}^{\left(\tau +1\right)}={\hat{b}}^{\left(\tau +1\right)}\left[a_0^{\left(\tau \right)}+{\delta }^{\left(\tau \right)}-P\right]+1\\ ={\hat{b}}^{\left(\tau +1\right)}\left[a_0^{\left(\tau +1\right)}-P\right]+1\end{array},& \left(18\right)\end{array}$

where a₀^(τ)=Σ_i^P a_i^(τ) and a₀^(τ)+δ^(τ). Roughly speaking, the larger the a₀^(τ), the smaller the total variance of the Dirichlet distribution. Adding δ^(τ)≥0 leads to a₀^(τ+1)≥a₀^(τ), such that the variance reduces monotonically with each iteration τ. Note that the introduction of δ^(τ) does not change the position of the mode in the distribution. We set δ^(τ) to a constant for all τ.

We demonstrate the performance of the proposed adaptive sequential measurement scheme for 10 distinct realizations of 4 point emitter constellations. The position of the point emitters are generated inside a circle with radius of 0.25 rl (Rayleigh length) with the minimum emitter separation set to 0.1 rl. The relative brightness are chosen to be uniform. The total photon budget is set to 5×10⁵ and each adaptive sequential step utilizes 10⁴ photons. The adaptive sequential scheme is initialized by employing 1000 photons for a direct imaging measurement followed by expectation maximization (EM) algorithm to estimate the initial model parameters. The remaining photons are detected by using Personick projection measurement in each adaptive sequential step.

For traditional direct imaging (baseline), Richardson-Lucy deconvolution algorithm is first used deconvolve the blurred image followed by the k-mean clustering algorithm to find the position and relative brightness of identified point emitters. If P is known a priori, the k-mean clustering algorithm is used once for P; otherwise the algorithm is used for all P_min≤P≤P_max and the model (n) that produces the smallest average (over all P) position error is picked. Here the position error for the i^th point emitter is defined as:

$\sqrt{{\left(x_i-\widehat{x_i}\right)}^2+{\left(y_i-\widehat{y_i}\right)}^2}.$

For each of the 10 constellations, we employ 100 Monte Carlo simulation. FIG. 5A shows an illustrative realization of the point emitter cluster and estimated location and brightness using the two measurement schemes. The example of FIG. 5A corresponds to estimation given by direct imaging and Personick projection measurement, when the number of emitters is given as a prior. The solid dots, circles, and squares are the ground truth, estimation from direct imaging, and Personick projection measurement respectively. All the marker sizes are proportional to the brightness. To obtain the average performance of the proposed measurement scheme, for each point emitter realization, we first pair the ground truth point emitter location with the estimated locations, such that the sum of the position errors over all point-source matched pairs is minimized. The average (over all emitters) position error distribution of the point emitters is shown in FIG. 5B. We observe that the proposed adaptive scheme outperforms the direct imaging. More specifically, the mean position error obtained by the adaptive scheme is four-fold lower than that of the direct imaging.

When P_max=6 is given as a prior, the estimation algorithm has to also estimate P. For the same set of constellations and using the same number of simulations, the distribution of number of point, emitters estimated by the two measurement schemes in shown in FIG. 5C. We observe that the adaptive Personick projective scheme estimates the correct number of point emitters with a 30% success rate relative to only 10% for direct imaging. FIG. 5D shows the corresponding position error distribution, computed only for cases where P is estimated correctly. We observe that when the P is estimated correctly the proposed adaptive scheme maintains the significant performance advantage over direct imaging.

Based on quantum estimation theory, we have described examples of a quantum adaptive Bayesian multi-parameter estimation scheme. In prophetic examples, we applied our proposed approach to super-resolving point emitters in an imaging application and demonstrated its superiority to the state-of the-art direct imaging approach. This demonstrates that quantum estimation theory applied to sensing and imaging measurement design problems can yield significant gains by increasing (in some cases maximizing) the photon information extraction with certain optical measurements.

The spatial mode sorting may be performed with various optical configurations, as discussed below.

FIG. 6A shows an example of a spatial anode sorting system 600. The incoming beam 602 reflects off a spatial light modulator 604 containing five independently controlled regions that modifies the intensity or phase of the incoming beam 602. A mirror 606 reflects the incoming beam 602 after it has interacted with one or more of the regions of the spatial light modulator 604. The incoming beam 602 is then sent to a detector 608, such as an EMCCD or CMOS camera. The information produced by the detector is then sent to a processor 610, such as an FPGA, which can then control the intensity and phase of future incoming light 602 after reflecting from the spatial light modulator 604.

FIG. 6B shows the spatial mode sorting system 600 with an incoming beam containing a first mode 620A and sorting it to a first region of the detector image 622A.

FIG. 6C shows the spatial anode sorting system 600 with an incoming beam containing a second mode 620B and sorting it to a second region of the detector image 622B.

FIG. 6D shows the spatial mode sorting system 600 with an incoming beam containing a first mode 620C and sorting it to a third region of the detector image 622C.

If a superposition of the modes 620A, 620B, and 620C is received in the beam 602, the ratio of the spot intensities on the resulting detector image can be used to infer the relative strength of the modes in the received beam 602.

FIG. 7 shows a second example of a spatial mode sorting system 700. A first spatial light modulator 701 reflects and modifies the intensity or phase of the incoming beam 710. A second spatial light modulator 702, a third spatial light modulator 703, a fourth spatial light modulator 704, and a fifth spatial light modulator 705 further reflect and modify the intensity or phase of the incoming beam 710.

FIG. 8 shows a third example of a spatial mode sorting system 800. A first spatial light modulator 801 transmits and modifies the intensity or phase of the incoming beam 810. A second spatial light modulator 802, a third spatial light modulator 803, a fourth spatial light rrrodulator 804, and a fifth spatial light modulator 805 further transmit and modify the intensity or phase of the incoming beam 810.

The techniques described above for controlling and configuring a spatial mode sorting system can be implemented using software for execution on a computer system. For example, the software can define procedures in one or more computer programs that execute on one or more programmed or programmable computer systems (e.g., desktop, distributed, client/server computer systems) each including at least one processor, at least one data storage system (e.g., including volatile and non-volatile memory and/or storage elements), at least one input device (e.g., keyboard and mouse) or port, and at least one output device (e.g., monitor) or port. The software may form one or more modules of a larger program.

The software may be provided on a non-transitory medium such as a computer-readable storage medium (e.g., solid state memory or media, or magnetic or optical media) readable by a general or special purpose programmable computer system, or delivered over a communication medium (e.g., encoded in a propagated signal) such as network to a computer system where it is stored in a non-transitory medium and executed. Each such computer program can be used to configure and operate the computer system when the non-transitory medium is read by the computer system to perform the procedures of the software.

While the disclosure has been described in connection with certain embodiments, it is to be understood that the disclosure is not to be limited to the disclosed embodiments but, on the contrary, is intended to cover various modifications and equivalent arrangements included within the scope of the appended claims, which scope is to be accorded the broadest interpretation so as to encompass all such modifications and equivalent structures as is permitted under the law.

Claims

1. A method for imaging a distribution of one or more optical sources, the method comprising: receiving respective optical signals from a spatial mode sorter during each of two or more detection intervals of time;after each of the two or more detection intervals of time, processing information based at least in part on: (1) the respective optical signal received in the corresponding detection interval of time, and (2) a first set of models comprising a set of distributions related to one or more optical sources, each model corresponding to a different number of optical sources in the distribution, andconfiguring the spatial mode sorter based at least in part on the processing; andproviding an estimated measurement characterizing the distribution of one or more optical sources based at least in part on the processed information;wherein the processing after at least one of the two or more detection intervals of time includes computing an eigen-projection.

2. The method of claim 1, wherein the configuring after each of the two or more detection intervals of time configures the spatial mode sorter to project the respective optical signals onto a basis of the computed eigen-projection.

3. The method of claim 1, wherein the first set of models includes a set of spatial distributions for the optical sources.

4. The method of claim 3, wherein the first set of models further includes a set of Bayesian prior probability distributions for the set of spatial distributions.

5. The method of claim 4, wherein the set of Bayesian prior probability distributions for the set of spatial distributions includes a set of Gaussian distributions.

6. The method of claim 5, wherein a set of hyper-parameters for the Gaussian distributions are based at least in part on a result of an expectation maximization calculation that is based at least in part on the processed information.

7. The method of claim 1, wherein the first set of models includes a set of brightness distributions for the optical sources.

8. The method of claim 7, wherein the first set of models further includes a set of Bayesian prior probability distributions for the set of brightness distributions.

9. The method of claim 8, wherein the set of Bayesian prior probability distributions for the set of brightness distributions includes a set of Dirichlet distributions.

10. The method of claim 9, wherein a set of hyper-parameters for the Dirichlet distribution is based at least in part on a result of expectation maximization calculation that is based at least in part on the processed information.

11. The method of claim 1, wherein the eigen-projection is the eigenvectors of a symmetric logarithmic derivative operator.

12. The method of claim 11, wherein the symmetric logarithmic derivative operator is based at least in part on a set of one or more operators constructed from a single-parameter inference setting.

13. The method of claim 1, wherein the eigen-projection is the eigenvectors of an operator constructed from a Bayesian inference setting.

14. The method of claim 13, wherein the operator constructed from a Bayesian inference setting is based at least in part on a set of one or more operators constructed from a single-parameter Bayesian inference setting.

15. The method of claim 1, wherein the eigen-projection comprises a Personick eigen-projection.

16. The method of claim 1, wherein the processed information includes a second set of models.

17. The method of claim 16, the second set of models comprising a second set of distributions related) one or more optical sources, determined by (1) the first set of distributions related to one or more optical sources and (2) respective optical signals received in a previous detection interval of time.

18. The method of claim 1, wherein the processing after at least one of the two or more detection intervals of time includes computing a quantum Fisher information matrix associated with the respective optical signals.

19. The method of claim 1, wherein the processing after at least one of the two or more detection intervals of time includes computing a modified quantum Fisher information matrix derived in a Bayesian inference setting and associated with the respective optical signals.

20. One or more non-transitory computer-readable media, having instructions stored thereon that, when executed by a computer system, cause the computer system to perform operations comprising: receiving respective optical signals from a spatial mode sorter during each of two or more detection intervals of time;after each of the two or more detection intervals of time, processing information based at least in part on: (1) the respective optical signal received in the corresponding detection interval of time, and (2) a first set of models comprising a set of distributions related to one or more optical sources, each model corresponding to a different number of optical sources in the distribution, andconfiguring the spatial mode sorter based at least in part on the processing; andproviding an estimated measurement characterizing the distribution of one or more optical sources based at least in part on the processed information;wherein the processing after at least one of the two or more detection intervals of time includes computing an eigen-projection.

21. An apparatus for imaging a distribution of one or more optical sources, the apparatus comprising: a spatial mode sorter that defines a configurable basis comprising a set of spatial modes onto which an incoming optical signal is projected; anda control module configured to: receive respective optical signals from the spatial mode sorter during each of two or more detection intervals of time;after each of the two or more detection intervals of time, process information based at least in part on: (1) the respective optical signal received in the corresponding detection interval of time, and (2) a first set of models comprising a set of distributions related to one or more optical sources, each model corresponding to a different number of optical sources in the distribution, andconfigure the spatial mode sorter based at least in part on the processing; andprovide an estimated measurement characterizing the distribution of one or more optical sources based at least in part on the processed information;wherein he processing after at least one of the two or more detection intervals of time includes computing an eigen-projection.