arXiv CSopen access

Training Neural Networks to Approach the Optimum Bayes Estimator in Dense Multi-Emitter Localization

· arxiv_cs

arXiv CS · Papers · License: Open Access

machine learning, deep learning, neural networks

Training Neural Networks to Approach the Optimum Bayes Estimator in Dense MultiEmitter Localization YI SUN*, MONA SHARIFI, AND MUZNA YUMMAN Electrical Engineering Department Applied and Embodied AI Lab The City College of City University of New York New York, NY 10031, USA *Email: [email protected] ORCID: 0000-0002-2527-8311

Abstract: Multi-emitter localization from a single frame sets the density limit of single molecule localization microscopy (SMLM), i.e., the accuracy falls as the emitters overlap, whereas at severe overlap the Cramér-Rao bound diverges and the best unbiased estimator loses accuracy without limit. The optimum estimator in the mean square sense is the Bayes estimator that minimizes the risk over a prior of emitter configurations, whereas for two or more emitters it has no closed form and its evaluation needs the point spread function, the emitter intensity and the per-pixel noise, so a direct computation requires the full system model. We train neural networks on synthesized frames to approach this estimator. Five networks that span learnable activations, Kolmogorov-Arnold edges, convolution, and self attention are trained with a permutation-matched loss whose Bayes estimator is exactly the reported metric, then compared against the unbiased Gaussian information-achieving estimator for a frame (UGIA-F) and the expectation-maximization global maximum likelihood estimator (EM-GML). The networks estimate from a frame alone, whereas EM-GML uses the system parameters and UGIA-F uses the system parameters and the true positions, so among the three only UGIA-F is an oracle and the networks and EM-GML are practically useful. Simulation results show that on indistribution frames every network stays far below UGIA-F in mean square error where the emitters overlap and approaches EM-GML, e.g., the vision Transformer attains an average error of 31.6 nm against 30.3 nm for EM-GML and 172 nm for UGIA-F at a density of 56 emitters per µm2. The networks remain bounded where the oracle UGIA-F diverges, whereas they reach the accuracy of EM-GML, which confirms the hypothesis in a proof of concept on a field of view (FOV) of 700 × 700 nm2. The result justifies the future work on training neural networks to approach the optimum Bayes estimator on the large-size frames to achieve high-throughput large-FOV super spatiotemporal resolution SMLM. © 2026

Introduction

Single molecule localization microscopy (SMLM) [1-4] forms a super-resolution image by estimating the positions of individual fluorescent emitters from a sequence of diffraction limited frames. The resolution of the reconstruction and the temporal resolution, defined by the time needed to acquire it, are set against each other by the density of active emitters per frame, i.e., a sparse frame carries well separated emitters to be localized easily but many frames at a low temporal resolution are needed, whereas a dense frame shortens the acquisition, thus improving the temporal resolution, but places several overlapping emitters in a single point spread function (PSF). The estimation of severely overlapping emitters from one frame is therefore critical to achieve both super spatial and temporal resolutions, and it is the problem studied in this paper.

The accuracy of localization has long been analyzed through the Cramér-Rao bound (CRB), i.e., the lower bound on the covariance of any unbiased estimator, which for a single emitter gives the familiar dependence of the accuracy on the photon count and the background [5-7]. The accuracy depends on the PSF model, i.e., the Gaussian model in two dimensions [1] and the engineered functions that encode depth in three dimensions [9-12], and on the information content of the frame that sets the achievable segmentation and signal to noise ratio [13, 14]. The bound extends to several emitters through the Fisher information matrix of the multiemitter frame [15-17], and the unbiased Gaussian information achieving estimator for a frame (UGIA-F) realizes it as an oracle benchmark that attains the CRB on every frame [18, 19]. In counterpart, the global maximum likelihood (GML) estimator [19] is a practical estimator, which is approached by multi-start expectation maximization (EM) and is denoted EM-GML [19-21]. Both benchmark estimators have limits that matter exactly where the emitters overlap. The CRB diverges when two emitters merge, so the accuracy attainable under the unbiasedness constraint is unbounded and UGIA-F degrades without limit, whereas the GML estimator is asymptotically efficient yet possesses no optimality property at a finite photon count [19]. Both benchmarks also need to know the PSF, the emitter intensity and the per-pixel noise, and UGIAF needs the true positions in addition, so UGIA-F is an oracle that a real experiment cannot run whereas EM-GML uses only the system parameters and is practical. The optimum estimator in the mean square sense is neither of these. It is the Bayes estimator that minimizes the risk over a prior of emitter configurations [22], which is the estimator this paper takes as the target. The Bayes estimator carries a bias and is therefore free of the constraint that makes the CRB diverge, so it stays accurate where the unbiased benchmark does not. For two or more emitters, however, the matched form of the Bayes estimator has no closed form, and its evaluation needs the posterior and hence the full system model, so a direct computation again requires the full system model at every frame. A learned estimator, e.g., a neural network, that approximated the Bayes estimator would therefore be more accurate than the unbiased benchmark at overlap and, once trained, would evaluate the estimate in a single forward pass rather than by the per-frame posterior integration that a direct computation requires. Neural networks have been applied to SMLM as learned localizers that map a frame directly to emitter positions after training on simulated frames, e.g., the convolutional network DeepSTORM [23], the dense-emitter network DECODE [24], the emission-pattern network smNet [25], the three-dimensional network DeepSTORM3D [26], and the residual deconvolutional network [27]. They remove the per-frame likelihood optimization and run in a single forward pass and they reach high emitter density and high speed. However, their relation to the estimators of estimation theory has not been made explicit, i.e., they are trained and evaluated empirically and are not formulated within the Bayes estimator framework, so it is not established which estimator a frame-trained network approaches, whether the network is optimal in any defined sense, or how it compares with the benchmark estimators of estimation theory. This paper addresses these questions that are theoretically and practically critical when applying neural networks to emitter localization. The contribution of this paper is twofold. First, we show that a network trained by minimizing an average loss over frames drawn from a prior approaches the Bayes estimator of that loss, and that training with a permutation-matched loss makes the target the matched Bayes estimator whose risk is exactly the localization metric reported here, so the choice of the training loss determines the obtained estimator in the Bayse sense. Second, we test by simulation the hypothesis that such a network outperforms the unbiased oracle UGIA-F and approaches the likelihood estimator EM-GML. Five networks that span learnable activations, Kolmogorov-Arnold edges, convolution, and self attention are trained under a common budget, so that a result common to all of them is attributed to learning from frames rather than to any single design. The simulation confirms the hypothesis on in-distribution frames and locates the advantage of the networks in the overlapped configurations where the oracle UGIA-F diverges.

The networks are trained on frames synthesized with the same system model that EM-GML and UGIA-F use, so the three estimators are built with the same system knowledge, and the networks additionally use the prior over configurations; among the three only UGIA-F uses the true positions and is an oracle, whereas the networks and EM-GML are practically useful, so the accuracy of the networks reflects learning the Bayes-optimal estimator rather than any information advantage. The study is a proof of concept carried out on a field of view (FOV) of 700 × 700 nm2. The results and findings motivate and justify the future work to train neural networks on the practically large frames to achieve high-throughput large-FOV super spatiotemporal resolution SMLM. The paper is organized as follows. Section 2 defines the frame model and the estimators through the Bayes risk framework, i.e., the minimum mean square error (MMSE), the sorted Bayes, the matched Bayes, the GML and the UGIA-F estimators, together with the localization metric. Section 3 defines the neural network estimators and establishes which estimator a frame-trained network approaches. Section 4 computes the matched Bayes estimator by quadrature for two emitters, so that the networks can be measured against the optimum. Section 5 reports the simulation results on an in-distribution dataset and on out-of-distribution circle constellations, together with a stratified analysis by emitter overlap and a verification of the EM-GML benchmark. Section 6 discusses the findings and Section 7 concludes. 2.

Estimators and the Bayes risk

2.1 Data frame model and log-likelihood The data frame model of this paper is the universal model of Refs. [15, 17, 19] that is broadly applicable in practical systems and experiments. For simplicity, considered are twodimensional (2D) imaging, a known emitter number and a known emitter intensity that is the same for all emitters. The analysis and results can be straightforwardly extended to other settings. The model and the estimators are stated here in the notation of Ref. [19] so that the results of the two studies can be read together. Consider 𝑀𝑀 emitters activated in a data frame. The 𝑚𝑚th emitter is positioned at 𝜽𝜽𝑚𝑚 = (𝑥𝑥𝑚𝑚 , 𝑦𝑦𝑚𝑚 )T ∈ ℝ2 in nm and all positions are stacked in the 2𝑀𝑀 dimensional vector 𝜽𝜽 = (𝜽𝜽1T , … , 𝜽𝜽T𝑀𝑀 )T .

(1)

𝑠𝑠(𝒌𝒌) = Δ𝑡𝑡 Δ𝑥𝑥 Δ𝑦𝑦 𝐼𝐼 𝑄𝑄(𝒌𝒌)

(2)

Every emitter emits 𝐼𝐼 photons per second on average and 𝐼𝐼 is known, so that 𝛽𝛽𝑚𝑚 = 1 for all 𝑚𝑚 in the notation of Ref. [17, 19]. The camera has 𝐾𝐾𝑥𝑥 × 𝐾𝐾𝑦𝑦 pixels of sizes Δ𝑥𝑥 and Δ𝑦𝑦 in nm, i.e., the FOV is [0, 𝐿𝐿𝑥𝑥 ] × [0, 𝐿𝐿𝑦𝑦 ] = [0, 𝐾𝐾𝑥𝑥 Δ𝑥𝑥 ] × [0, 𝐾𝐾𝑦𝑦 Δ𝑦𝑦 ], the frame time is Δ𝑡𝑡 in second, the pixel index is 𝒌𝒌 = �𝑘𝑘𝑥𝑥 , 𝑘𝑘𝑦𝑦 � and the pixel set is Ω. A photon emitted from the 𝑚𝑚th emitter arrives in the camera plane with the PSF, i.e., the probability density funciton 𝑞𝑞𝑚𝑚 (𝒖𝒖), 𝒖𝒖 ∈ ℝ2 , whose average over the 𝒌𝒌th pixel is 𝑞𝑞𝑚𝑚 (𝒌𝒌), 𝒌𝒌 ∈ Ω. The signal mean, i.e., the mean photon count from all emitters, in the 𝒌𝒌th pixel is where

𝑀𝑀

𝑄𝑄(𝒌𝒌) = � 𝛽𝛽𝑚𝑚 𝑞𝑞𝑚𝑚 (𝒌𝒌).

(3)

𝑏𝑏(𝒌𝒌) = Δ𝑡𝑡 Δ𝑥𝑥 Δ𝑦𝑦 𝑏𝑏𝒌𝒌 .

(4)

𝑚𝑚=1

The background autofluorescence and the camera readout together produce the mixed noise whose spatiotemporal density in photons per nm2 per second is 𝑏𝑏𝒌𝒌 and whose mean photon count in the 𝒌𝒌th pixel is

The mean of the Gaussian readout does not affect the Fisher information and a Gaussian variable whose mean equals its variance is well approximated by a Poisson variable and conversely [15]. The readout mean is therefore set equal to its variance, so that the mixed noise is Poisson and 𝑏𝑏𝒌𝒌 denotes the combined background and readout density. Ref. [19] adopts this approximation so that the EM algorithm operates on a Poisson likelihood, and it is adopted here for the same reason. Unlike the earlier studies [15, 17] where 𝑏𝑏𝒌𝒌 is spatially uniform, here 𝑏𝑏𝒌𝒌 varies from pixel to pixel, which models the slowly spatially varying autofluorescence and the pixel dependent readout of an sCMOS camera [28]. The pixel value is the sum of the signal and the noise, with

𝑉𝑉(𝒌𝒌) = 𝑆𝑆(𝒌𝒌) + 𝐵𝐵(𝒌𝒌) ∼ Poisson�𝑣𝑣(𝒌𝒌)�

(5)

𝑣𝑣(𝒌𝒌) = 𝑠𝑠(𝒌𝒌) + 𝑏𝑏(𝒌𝒌)

(6)

where 𝑆𝑆(𝒌𝒌) and 𝐵𝐵(𝒌𝒌) are mutually independent and the pixel values 𝑉𝑉(𝒌𝒌) are independent across pixels. The data frame is denoted by 𝑉𝑉 and its probability mass function is 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) = � 𝒌𝒌∈Ω

𝑣𝑣(𝒌𝒌)𝑉𝑉(𝒌𝒌) 𝑒𝑒 −𝑣𝑣(𝒌𝒌) , 𝑉𝑉(𝒌𝒌)!

(7)

which is also the likelihood function of 𝜽𝜽. Dropping the constant that is independent of 𝜽𝜽, the per-frame log-likelihood is ℓ1 (𝜽𝜽) = �[𝑉𝑉(𝒌𝒌) ln 𝑣𝑣 (𝒌𝒌) − 𝑣𝑣(𝒌𝒌)]

(8)

𝒌𝒌∈Ω

and the Fisher information matrix [15, 19] is 𝐅𝐅(𝜽𝜽) = �

𝒌𝒌∈Ω

1 𝜕𝜕𝜕𝜕(𝒌𝒌) 𝜕𝜕𝜕𝜕(𝒌𝒌) . 𝑣𝑣(𝒌𝒌) 𝜕𝜕𝜽𝜽 𝜕𝜕𝜽𝜽T

(9)

An estimator estimates 𝜽𝜽 from the single frame 𝑉𝑉, that is, an estimator is a function 𝜽𝜽ˆ (𝑉𝑉). 2.2 Bayes risk

The estimators compared in this paper are constructed on different principles, so a common criterion is needed before any of them can be judged. The criterion adopted here is the Bayes risk, which is the average estimation error over both the noise in the frame and the emitter positions that produced it. Let 𝑝𝑝(𝜽𝜽) denote the prior, that is, the distribution from which the emitter positions are drawn. When 𝜽𝜽 is regarded as random with prior 𝑝𝑝(𝜽𝜽), the function 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) of Eq. (7) is the conditional probability mass function of 𝑉𝑉 given 𝜽𝜽, so the joint distribution of the frame and the emitter positions is 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽)𝑝𝑝(𝜽𝜽). Let 𝐿𝐿�𝜽𝜽ˆ , 𝜽𝜽� denote a loss function that measures the error of an estimate. The risk of an estimator is its loss averaged over the joint distribution, 𝑅𝑅�𝜽𝜽ˆ � = � � 𝐿𝐿 �𝜽𝜽ˆ (𝑉𝑉), 𝜽𝜽� 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) 𝑝𝑝(𝜽𝜽) 𝑑𝑑𝑑𝑑 𝑑𝑑𝜽𝜽.

(10)

𝑝𝑝(𝑉𝑉) = � 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) 𝑝𝑝(𝜽𝜽) 𝑑𝑑𝜽𝜽

(11)

The Bayes estimator is defined as the estimator that minimizes 𝑅𝑅 over all functions 𝜽𝜽ˆ (⋅) for the loss 𝐿𝐿 and the prior 𝑝𝑝. Its risk is the smallest attainable by any estimator whatsoever, so it is the bottom line against which every other estimator can be measured. The Bayes estimator has a simple characterization. Let

be the marginal distribution of the frame, that is, the distribution of a frame produced by an emitter configuration drawn from the prior. The posterior of the emitter positions given the observed frame is then 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉) =

so that the joint distribution factors as

𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) 𝑝𝑝(𝜽𝜽) , 𝑝𝑝(𝑉𝑉)

(12)

𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) 𝑝𝑝(𝜽𝜽) = 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉) 𝑝𝑝(𝑉𝑉).

(13)

𝑅𝑅�𝜽𝜽ˆ � = � �� 𝐿𝐿 �𝜽𝜽ˆ (𝑉𝑉), 𝜽𝜽� 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉) 𝑑𝑑𝜽𝜽� 𝑝𝑝(𝑉𝑉) 𝑑𝑑𝑑𝑑.

(14)

𝜽𝜽ˆ Bayes (𝑉𝑉) = arg min � 𝐿𝐿 (𝒖𝒖, 𝜽𝜽) 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉) 𝑑𝑑𝜽𝜽,

(15)

Substituting Eq. (13) into Eq. (10) and integrating over 𝜽𝜽 first gives

The inner integral of Eq. (14) depends on the estimator only through the value 𝜽𝜽ˆ (𝑉𝑉) that the estimator takes at the observed frame, and 𝑝𝑝(𝑉𝑉) is non-negative, so the risk is minimized by minimizing the inner integral separately for every 𝑉𝑉. The Bayes estimator is therefore 𝒖𝒖

which is computed frame by frame although it minimizes an average over frames. Two consequences are used throughout this paper. First, the Bayes estimator depends on the loss function, so different losses define different optimal estimators even for the same data. Second, the optimality holds only when the prior in the posterior is the distribution the frames were actually drawn from, so the guarantee is lost when an estimator built for one prior is applied to data generated by another. Sections 2.3 to 2.5 derive the Bayes estimator for three loss functions. Each of the three is a loss with which a neural network can be trained, so the three subsections together state which estimator a network approaches for a given choice of the training loss, and thereby the principle by which a network is to be trained. Sections 2.6 and 2.7 then present the GML and UGIA-F estimators. Neither of the two is a Bayes estimator, and both are used in this paper only as benchmarks against which the networks are compared. 2.3 Minimum mean square error estimator The natural loss for localization is the squared error 𝑀𝑀

𝐿𝐿SE (𝒖𝒖, 𝜽𝜽) = � ‖𝒖𝒖𝑚𝑚 − 𝜽𝜽𝑚𝑚 ‖2 ,

(16)

𝜽𝜽ˆ MMSE (𝑉𝑉) = 𝔼𝔼[𝜽𝜽 ∣ 𝑉𝑉].

(17)

𝑚𝑚=1

so the Bayes estimator of Eq. (15) minimizes the mean square error (MSE), that is 𝔼𝔼[𝐿𝐿SE (𝒖𝒖, 𝜽𝜽) ∣ 𝑉𝑉] = 𝔼𝔼[‖𝒖𝒖 − 𝜽𝜽‖2 ∣ 𝑉𝑉] , where the expectation is taken with respect to 𝜽𝜽 conditioned on 𝑉𝑉. It is easy to obtain that the minimum mean square error (MMSE) estimator is the conditional mean This estimator is degenerate for multiple emitters, as the following proposition shows.

Proposition 1: Let 𝒮𝒮𝑀𝑀 be the set of permutations of {1, … , 𝑀𝑀} and for 𝜋𝜋 ∈ 𝒮𝒮𝑀𝑀 let 𝐏𝐏𝜋𝜋 permute the emitter blocks of 𝜽𝜽. If the prior is exchangeable, that is 𝑝𝑝(𝐏𝐏𝜋𝜋 𝜽𝜽) = 𝑝𝑝(𝜽𝜽) for all 𝜋𝜋 ∈ 𝒮𝒮𝑀𝑀 , then 𝔼𝔼[𝜽𝜽1 ∣ 𝑉𝑉] = 𝔼𝔼[𝜽𝜽2 ∣ 𝑉𝑉] = ⋯ = 𝔼𝔼[𝜽𝜽𝑀𝑀 ∣ 𝑉𝑉].

(18)

is yielded by the MMSE estimator.

█

The proof is given in Appendix A.1. Since the emitters in an experiment are physically indistinguishable, the prior is exchangeable whenever the emitter positions are drawn independently, so Proposition 1 applies to the setting of this paper, i.e., emitter positions drawn independently, and to any other setting in which the emitters are not labeled by some further observable. All 𝑀𝑀 components of 𝜽𝜽ˆ MMSE coincide at the common posterior mean of Eq. (18), so the MMSE estimator reports a single point regardless of how the emitters are arranged and it can never resolve them. The reason for the degeneracy is the mismatch between what the loss requires and what the frame supplies. The loss 𝐿𝐿SE of Eq. (16) is defined with respect to the emitter indices, since the 𝑚𝑚th estimate is penalized by its distance to the 𝑚𝑚th emitter and to no other. The frame however carries no information that distinguishes one index from another, because the likelihood of Eq. (7) is invariant under any permutation of the indices. The posterior marginal of 𝜽𝜽𝑚𝑚 is therefore the same for every 𝑚𝑚, and the value of 𝒖𝒖𝑚𝑚 that minimizes 𝔼𝔼[‖𝒖𝒖𝑚𝑚 − 𝜽𝜽𝑚𝑚 ‖2 ∣ 𝑉𝑉] is that one common posterior mean for every 𝑚𝑚. In short, the loss demands a correspondence between estimates and emitters that the posterior does not provide, and the estimator answers by giving the same estimate to all of them. For 𝑀𝑀 = 2 the resulting MMSE has a closed form. Corollary 1: Let 𝑀𝑀 = 2. Denote by 𝑑𝑑 = ‖𝜽𝜽1 − 𝜽𝜽2 ‖ the emitter separation, 𝒄𝒄 = (𝜽𝜽1 + 𝜽𝜽2 )/ ˆ the common posterior mean of Eq. (18). Then 2 the midpoint of the two emitters and 𝒎𝒎 2

is the MMSE of Eq. (16)

ˆ − 𝜽𝜽𝑚𝑚 ‖2 = � ‖𝒎𝒎

𝑚𝑚=1

𝑑𝑑 2 ˆ − 𝒄𝒄‖2 + 2‖𝒎𝒎 2

(19)

█

The proof is given in Appendix A.2. The degenerate estimator is therefore the centroid estimator, whose error is set by half the emitter separation and is small only when the emitters are unresolvable. 2.4 The sorted Bayes estimator The degeneracy of Section 2.3 follows from the exchangeability of the posterior, so it can be removed by breaking that exchangeability. The simplest way is to impose a canonical order on the emitters, for instance by their 𝑥𝑥 coordinates. Let 𝜽𝜽(1) , … , 𝜽𝜽(𝑀𝑀) denote the emitter positions arranged so that 𝑥𝑥(1) ≤ ⋯ ≤ 𝑥𝑥(𝑀𝑀) and consider the sorted squared error 𝑀𝑀

𝐿𝐿S (𝒖𝒖, 𝜽𝜽) = � �𝒖𝒖𝑚𝑚 − 𝜽𝜽(𝑚𝑚) � .

(20)

𝜽𝜽ˆ S (𝑉𝑉) = 𝔼𝔼�𝜽𝜽(⋅) ∣ 𝑉𝑉�,

(21)

𝑚𝑚=1

The prior of the sorted vector is not exchangeable, so Proposition 1 no longer applies. It is easy to obtain that the Bayes estimator of Eq. (15) for this loss is returning 𝑀𝑀 distinct estimates. It is the posterior mean of the order statistics of the emitter positions along the 𝑥𝑥 axis. The sorting removes the degeneracy at a price, which is an anisotropy between the sorting coordinate 𝑥𝑥 and the coordinate 𝑦𝑦 orthogonal to it. The order is decided by the 𝑥𝑥 coordinates alone, so whenever two emitters have close 𝑥𝑥 coordinates, the noise decides which of them is called the first, and the 𝑦𝑦 estimate that accompanies the decision jumps between the two emitters. The following proposition quantifies the effect for two emitters.

Proposition 2: Let 𝑀𝑀 = 2 and let the two estimates be ordered by their 𝑥𝑥 coordinates. For a fixed true configuration let the unordered estimates over repeated frames be 𝑋𝑋𝑚𝑚 = 𝑥𝑥𝑚𝑚 + 𝜉𝜉𝑚𝑚 and 𝑌𝑌𝑚𝑚 = 𝑦𝑦𝑚𝑚 + 𝜂𝜂𝑚𝑚 for 𝑚𝑚 = 1, 2, where 𝜉𝜉1 , 𝜉𝜉2 , 𝜂𝜂1 , 𝜂𝜂2 are mutually independent and 𝒩𝒩(0, 𝜎𝜎 2 ). Write 𝑑𝑑𝑥𝑥 = |𝑥𝑥1 − 𝑥𝑥2 | and 𝑑𝑑𝑦𝑦 = |𝑦𝑦1 − 𝑦𝑦2 | for the coordinate separations of the two emitters and thus 𝑝𝑝 = Φ �−

𝑑𝑑𝑥𝑥

𝜎𝜎√2

(22)

�

is the probability that the order of 𝑋𝑋1 and 𝑋𝑋2 disagrees with the order of 𝑥𝑥1 and 𝑥𝑥2 , where Φ is the standard normal distribution function. Then the ordered estimate of the 𝑦𝑦 coordinate has variance

(23)

Var�𝑌𝑌(1) � = 𝜎𝜎 2 + 𝑝𝑝(1 − 𝑝𝑝) 𝑑𝑑𝑦𝑦2

whereas the ordered estimate of the 𝑥𝑥 coordinate has variance determined by the mean absolute difference of the two estimates,

namely

𝔼𝔼 [|𝑋𝑋1 − 𝑋𝑋2 |] =

2𝜎𝜎

√𝜋𝜋

exp �−

𝑑𝑑𝑥𝑥2 � + 𝑑𝑑𝑥𝑥 (1 − 2𝑝𝑝), 4𝜎𝜎 2

1 Var�𝑋𝑋(1) � = 𝜎𝜎 2 + [𝑑𝑑𝑥𝑥2 − 𝔼𝔼2 (|𝑋𝑋1 − 𝑋𝑋2 |)], 4

which increases monotonically from 𝜎𝜎 2 (1 − 1/𝜋𝜋) at 𝑑𝑑𝑥𝑥 = 0 to 𝜎𝜎 2 as 𝑑𝑑𝑥𝑥 → ∞.

(24)

(25) █

The proof is given in Appendix A.3. Two consequences follow. The excess variance 𝑝𝑝(1 − 𝑝𝑝)𝑑𝑑𝑦𝑦2 appears only in the coordinate that is not used for sorting, so the estimates of one emitter obtained from repeated frames scatter more widely along 𝑦𝑦 than along 𝑥𝑥. And the excess is largest when the two emitters are aligned perpendicular to the sorting axis, i.e., 𝑑𝑑𝑥𝑥 ≪ 𝑑𝑑𝑦𝑦 . Since 𝑑𝑑𝑥𝑥 = 0 gives 𝑝𝑝 = 1/2 and Var�𝑌𝑌(1) � = 𝜎𝜎 2 +

𝑑𝑑𝑦𝑦2 , 4

(26)

the 𝑦𝑦 error, i.e., the square root of the variance of Eq. (26), approaches 𝑑𝑑𝑦𝑦 /2, which is the error of the degenerate MMSE estimator in Corollary 1. In that geometry sorting therefore recovers nothing beyond the degenerate collapse that it was introduced to avoid. Conversely, when 𝑑𝑑𝑥𝑥 is large compared with 𝜎𝜎, the probability 𝑝𝑝 vanishes and the sorted estimator agrees with the matched estimator of Section 2.5, implying the sorting on the 𝑥𝑥 coordinates alone is sufficient. For 𝑀𝑀 > 2 the same anisotropy arises for every pair of emitters whose 𝑥𝑥 coordinates are close. Therefore, sorting on the 𝑥𝑥 coordinates alone still causes ambiguity in the estimates of 𝑦𝑦 coordinates particularly when the 𝑥𝑥 coordinates are close. 2.5 The matched Bayes estimator

The sorted loss of Eq. (20) breaks the exchangeability by an arbitrary convention rather than by the geometry of the estimate, which is the origin of the anisotropy of Proposition 2. The alternative is to let the correspondence between estimates and emitters be chosen so as to minimize the error itself, which gives the permutation matched squared error

whose Bayes estimator is

𝑀𝑀

𝐿𝐿M (𝒖𝒖, 𝜽𝜽) = min � �𝒖𝒖𝑚𝑚 − 𝜽𝜽𝜋𝜋(𝑚𝑚) � , 𝜋𝜋∈𝒮𝒮𝑀𝑀

𝑚𝑚=1

(27)

𝑀𝑀

2 𝜽𝜽ˆ B (𝑉𝑉) = arg min � � min � �𝒖𝒖𝑚𝑚 − 𝜽𝜽𝜋𝜋(𝑚𝑚) � � 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉) 𝑑𝑑𝜽𝜽. 𝒖𝒖

𝜋𝜋∈𝒮𝒮𝑀𝑀

(28)

𝑚𝑚=1

Here 𝜋𝜋(𝑚𝑚) is the index of the true emitter matched to the 𝑚𝑚th estimate by the permutation 𝜋𝜋 ∈ 𝒮𝒮𝑀𝑀 . The minimization over 𝜋𝜋 inside the integral removes the dependence of the loss on the emitter indices, so the exchangeability of the posterior no longer forces the estimates together, and 𝜽𝜽ˆ B returns 𝑀𝑀 distinct estimates. The minimizing permutation is the Hungarian assignment [29]. This is the reason why the Hungarian assignment is required in the training of a network. Without it the training loss is 𝐿𝐿SE of Eq. (16), whose Bayes estimator is the degenerate single posterior mean of Proposition 1, so a network trained without the assignment converges to 𝑀𝑀 coincident estimates and can never resolve the emitters. The sorting on the 𝑥𝑥 coordinates alone also yields large variances on the estimates of 𝑦𝑦 coordinates as indicated in Proposition 2. The Hungarian assignment is the sorting on both 𝑥𝑥 and 𝑦𝑦 coordinates, thus completely eliminating the index ambiguity. For 𝑀𝑀 = 1 the assignment is trivial and 𝜽𝜽ˆ B reduces to the posterior mean. 𝜽𝜽ˆ B has no closed form for 𝑀𝑀 ≥ 2 and is computed in Section 4 by quadrature over the posterior for 𝑀𝑀 = 2. It is the central object of this paper, since it is the minimizer of the risk under the loss that the quality metric of Section 2.8 measures, so no estimator can achieve a smaller value of that metric on frames drawn from the prior. In particular, it is not worse than either of the two estimators of Sections 2.3 and 2.4. The absence of a closed form is due to the minimization over permutations in Eq. (27), i.e., the winning permutation depends on the estimate so the objective of Eq. (28) is a non-convex piecewise quadratic function of the estimate whose minimizer cannot be written in closed form for 𝑀𝑀 ≥ 2. The estimate is instead obtained by the fixed-point iteration of Section 4, i.e., an assignment step that matches the posterior mass to the current estimate and an update step that sets each estimate to the posterior mean of the mass assigned to it, alternated to convergence. The iteration also needs the posterior 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉) of Eq. (12) and hence the likelihood of Eq. (7) with its PSF, its intensity and its per-pixel noise, so a direct computation of 𝜽𝜽ˆ B uses the full system model at every frame. This per-frame computational cost is one motivation for approximating 𝜽𝜽ˆ B by a network that, once trained, evaluates the estimate in a single forward pass, which is developed in Section 3. Corollary 2: Let 𝜽𝜽ˆ MMSE , 𝜽𝜽ˆ S and 𝜽𝜽ˆ B be the estimators of Eqs. (17), (21) and (28). Then 𝑅𝑅�𝜽𝜽ˆ B � ≤ 𝑅𝑅�𝜽𝜽ˆ MMSE �,

𝑅𝑅�𝜽𝜽ˆ B � ≤ 𝑅𝑅�𝜽𝜽ˆ S �

under the matched loss 𝐿𝐿M of Eq. (27).

(29)

█

The proof is given in Appendix A.4. The inequalities of Eq. (29) are not strict. By Proposition 2 the sorted estimator agrees with the matched estimator when the 𝑥𝑥 coordinates of every pair of emitters are well separated compared with 𝜎𝜎, so the two risks are equal in that limit and the gap opens only as emitters approach one another in the sorting coordinate. The corollary compares the risks, that is the averages over the prior, and it does not assert that 𝜽𝜽ˆ B is more accurate than the other two on every individual frame. 2.6 Global maximum likelihood estimator

The GML estimator [19] maximizes the per-frame log-likelihood, 𝜽𝜽ˆ GML = arg max ℓ1 (𝜽𝜽) . 𝜽𝜽

(30)

It is asymptotically consistent, asymptotically normal and asymptotically efficient as the Fisher information grows without bound [19], whereas it possesses no optimality property at a finite photon count. In the Bayes framework of Section 2.2, 𝜽𝜽ˆ GML is the mode of the posterior under a flat prior, that is under 𝑝𝑝(𝜽𝜽) = const, which is a different functional of the posterior from the conditional mean of Eq. (17). The expectation maximization (EM) algorithm [19] performs a local search that increases ℓ1 (𝜽𝜽) at every iteration and converges to a local maximum likelihood (LML) point. The likelihood function of multiple emitters may present several LML points, one of which is the global maximizer 𝜽𝜽ˆ GML , so a single EM run started at an arbitrary point converges to an LML point rather than to 𝜽𝜽ˆ GML . When the emitters overlap severely, the Fisher information matrix is near singular, so the log-likelihood is nearly flat about its maximizer and the ascent is slow, which further reduces the chance that a single run arrives at the global maximizer within a finite iteration budget. The global maximizer is therefore approached by multi-start EM. Let 𝜽𝜽ˆ (1) , … , 𝜽𝜽ˆ (𝐾𝐾) be the LML points to which 𝐾𝐾 EM runs converge from 𝐾𝐾 independent random initial positions and let 𝑘𝑘 ⋆ = arg max ℓ1 �𝜽𝜽ˆ (𝑘𝑘) �.

(31)

⋆ 𝜽𝜽ˆ EM-GML = 𝜽𝜽ˆ (𝑘𝑘 ) ,

(32)

1≤𝑘𝑘≤𝐾𝐾

The EM-GML estimate is the one that achieves the largest log-likelihood among the 𝐾𝐾 points,

which is taken as an approximation of 𝜽𝜽ˆ GML . Selection by Eq. (31) uses the log-likelihood alone, which depends only on the estimated positions and the observed frame, so the selection rule requires no knowledge of the true positions and a larger 𝐾𝐾 can only bring the estimate closer to 𝜽𝜽ˆ GML . The initial positions are drawn in a neighborhood of the true positions so that the multi-start search reliably reaches the global maximizer for the benchmark, whereas the global maximizer itself is a function of the frame and the system parameters and does not use the true positions, so EM-GML is a practical estimator. Section 5.6 verifies against an exhaustive search over the whole parameter space that 𝐾𝐾 = 250 attains the global maximizer on 99.7 percent of frames for 𝑀𝑀 = 2. In practice the true positions are unavailable and EM can be initialized by a successive interference cancellation (SIC) estimator [21] or any other estimator at the cost of losing the guarantee. 2.7 Unbiased Gaussian information achieving estimator The Cramér-Rao bound (CRB) is the matrix 𝐅𝐅(𝜽𝜽)−1 and the covariance of any unbiased estimator is bounded below by it in the Loewner ordering. The UGIA-F estimator [18, 19] is the Gaussian estimator 𝜽𝜽ˆ F ∼ 𝒩𝒩(𝜽𝜽, 𝐅𝐅(𝜽𝜽)−1 ),

(33)

𝜽𝜽ˆ F = 𝜽𝜽 + 𝐅𝐅(𝜽𝜽)−1/2 𝒈𝒈

(34)

which is unbiased and attains the Fisher information and the CRB. It is realized by where 𝐅𝐅(𝜽𝜽)−1/2 is the inverse square root of the Fisher information matrix and 𝒈𝒈 is a standard Gaussian vector. The UGIA-F estimator is an oracle benchmark rather than a deployable algorithm since its construction uses the true positions, and it represents the accuracy attainable by the best unbiased estimator. Unbiasedness is a constraint, not an optimality property. When two emitters approach each other, the Fisher information matrix becomes near singular and the CRB diverges, so the accuracy attainable under the unbiasedness constraint degrades without limit while a biased estimator such as the above EM-GML [19] remains free of that constraint. 2.8 Quality metric

Every estimator of Sections 2.3 to 2.7 produces 𝑀𝑀 estimated positions from one data frame, whereas the indices of the estimates are unknown, so the root mean square error (RMSE), 1/2 defined as �𝐿𝐿SE �𝜽𝜽ˆ , 𝜽𝜽�/𝑀𝑀� where 𝐿𝐿SE is the squared error of Eq. (16) evaluated under the true correspondence, cannot be evaluated and a correspondence between the estimated and the true positions must be established before an error can be computed. The correspondence is chosen by maximum likelihood. Under the model in which every estimate is Gaussian about its true position with a common standard deviation 𝜎𝜎, which holds approximately for any reasonable estimator [30], the likelihood of the correspondence 𝜋𝜋 is 𝑀𝑀

1 2 ln 𝐿𝐿 (𝜋𝜋) = const − 2 � �𝜽𝜽ˆ 𝜋𝜋(𝑚𝑚) − 𝜽𝜽𝑚𝑚 � , 2𝜎𝜎

(35)

𝑚𝑚=1

so the correspondence of maximum likelihood is the one that minimizes the total squared distance, which is the assignment computed by the Hungarian algorithm. The quality metric is therefore the root mean square minimum distance under that assignment, denoted RMSMD-H, 1/2

𝑀𝑀

1 2 𝐷𝐷H �𝜽𝜽ˆ , 𝜽𝜽� = � min � �𝜽𝜽ˆ 𝜋𝜋(𝑚𝑚) − 𝜽𝜽𝑚𝑚 � � 𝑀𝑀 𝜋𝜋∈𝒮𝒮𝑀𝑀 𝑚𝑚=1

(36)

Two properties of Eq. (36) are used repeatedly. First, 𝑀𝑀𝐷𝐷H2 is exactly the loss 𝐿𝐿M of Eq. (27), so the matched Bayes estimator 𝜽𝜽ˆ B of Eq. (28) is by construction the minimizer of the reported metric Eq. (36), and the training loss of Section 3.3 evaluates the same quantity. Second, since the minimization is over assignments, 𝐷𝐷H is not greater than the RMSE evaluated under the true correspondence where 𝜋𝜋 is an identity permutation. A data movie of 𝑁𝑁 frames is processed frame by frame, and the metric reported for the movie is the average RMSMD-H, denoted ARMSMD-H, 𝑁𝑁

1/2

1 𝐷𝐷H = � � 𝐷𝐷H2 �𝜽𝜽ˆ 𝑛𝑛 , 𝜽𝜽𝑛𝑛 �� 𝑁𝑁 𝑛𝑛=1

(37)

which is the Monte Carlo estimate of 𝑅𝑅�𝜽𝜽ˆ �/𝑀𝑀 under the loss 𝐿𝐿M . The averaging is performed on the squared distances and the square root is taken once at the end. The accuracy of an estimator reported in this paper is its RMSMD-H and ARMSMD-H of Eqs. (36) and (37), i.e., the smaller the value the more accurate. RMSMD-H and ARMSMD-H measure the average distance between the estimate and its true position on a 2D plane. When considering the average error in one axis as usually considered in microscopy, the metrics shall be divided by √2. The metric is defined frame by frame because this paper estimates the emitter positions from a single data frame, so every frame is an independent estimation problem in which the number of estimates equals the number of emitters. The root mean square minimum distance RMSMD [18] and its partition form RMSMD-P [30] evaluate a reconstructed SMLM image instead, in which the estimates of all frames of a movie are pooled and one emitter carries many estimates, so the number of estimates exceeds the number of emitters and no correspondence between the two sets exists. They are not used here. The companion study [19] reports the RMSE because there the indices of the estimates of both the UGIA-F estimator and the EMGML estimator are known, whereas here the neural networks produce unindexed estimates so the RMSE is unavailable and 𝐷𝐷H is used for all estimators so that they are compared on the same footing. 3.

Neural network estimators

3.1 What a frame-trained network estimates

A neural network estimator is a function 𝜽𝜽ˆ NN (𝑉𝑉; 𝒘𝒘) of the frame, parameterized by the network weights 𝒘𝒘. The weights are chosen to minimize the empirical average of a loss 𝐿𝐿 over a set of training frames, 𝑁𝑁t

1 𝒘𝒘 ˆ = arg min � 𝐿𝐿 �𝜽𝜽ˆ NN (𝑉𝑉𝑛𝑛 ; 𝒘𝒘), 𝜽𝜽𝑛𝑛 �, 𝒘𝒘 𝑁𝑁t

(38)

𝑛𝑛=1

where the 𝑁𝑁t training pairs (𝜽𝜽𝑛𝑛 , 𝑉𝑉𝑛𝑛 ) are generated by drawing 𝜽𝜽𝑛𝑛 from the prior 𝑝𝑝(𝜽𝜽) and then 𝑉𝑉𝑛𝑛 from the frame model of Eq. (5). Eq. (38) is a sample version of the risk 𝑅𝑅 of Eq. (10) under the loss 𝐿𝐿. As 𝑁𝑁t grows, the sample average converges to 𝑅𝑅, and the minimizer of 𝑅𝑅 over all functions of the frame is the Bayes estimator for that loss by Eq. (15). A frame-trained network therefore approximates the Bayes estimator of whichever loss it is trained with, and it departs from the Bayes estimator only through the finite capacity of the network, the finite number of training frames and the optimization. This is the theoretical basis on which a neural network is used as an estimator in this paper, and it is what distinguishes a network from the GML and UGIA-F estimators, which target the posterior mode and the best unbiased accuracy respectively. The three losses of Sections 2.3 to 2.5 therefore give three different networks from the same architecture and the same training frames. First, a network trained with the squared error 𝐿𝐿SE of Eq. (16) approaches the MMSE estimator of Eq. (17), which by Proposition 1 is degenerate, so its 𝑀𝑀 outputs converge to a single point and the emitters are never resolved. This holds whenever the emitter labels supplied in training are exchangeable, which is the case when the emitter positions are drawn independently and presented in the order generated. Second, a network trained with the sorted squared error 𝐿𝐿S of Eq. (20) approaches the sorted Bayes estimator of Eq. (21). Its outputs are distinct, whereas by Proposition 2 its estimates scatter more widely along the coordinate that is not used for sorting, and the excess is largest for emitter pairs aligned perpendicular to the sorting axis. Third, a network trained with the matched squared error 𝐿𝐿M of Eq. (27) approaches the matched Bayes estimator of Eq. (28). Its outputs are distinct and, by Corollary 2, its risk under the reported metric is not greater than that of either of the other two. The choice of the training loss therefore determines the estimator that is obtained, and the matched loss is chosen in this paper because it is the loss whose Bayes estimator minimizes exactly the metric of Section 2.8. The three networks are compared numerically in Section 5. Three further properties of Eq. (38) are worth stating explicitly. First, the network is built with the system parameters through its training frames. The training pairs of Eq. (38) are generated from the frame model of Eq. (5) with the PSF, the emitter intensity and the per-pixel noise, and from the true positions drawn from the prior, so the network is given the same system knowledge as the GML and UGIA-F estimators, i.e., it acquires that knowledge in training rather than by an explicit per-frame computation. Once trained, it maps a frame to positions in a single forward pass and does not recompute the likelihood. Second, the network additionally uses the prior. Training draws 𝜽𝜽𝑛𝑛 from 𝑝𝑝(𝜽𝜽), so the network learns the distribution of emitter configurations, which the GML and UGIA-F estimators do not use. The network therefore has the same system knowledge as EM-GML and UGIA-F and the prior in addition, which is the side information that the Bayes framework prescribes and is examined in Section 6. Third, the prior used in training is uniform. When the distribution of the emitter positions is unknown, the least informative choice on a bounded support ℬ is the uniform density, i.e., 𝑝𝑝(𝜽𝜽) constant on ℬ and zero outside, which is the maximum entropy prior on ℬ and commits to no particular arrangement of the emitters. Training on the uniform prior therefore makes the network accurate across the whole support rather than tuned to an assumed clustering, whereas a more informative prior can replace it whenever one is known.

Under the uniform prior, the estimator reduces to a likelihood weighted average. With 𝑝𝑝(𝜽𝜽) constant on ℬ, the posterior of Eq. (12) becomes the likelihood of Eq. (7) normalized over ℬ, i.e., for 𝜽𝜽 ∈ ℬ, 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉) =

𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) , ∫ℬ 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽′) 𝑑𝑑𝜽𝜽′

and the posterior is zero outside ℬ. The matched Bayes estimator of Eq. (28) is therefore 𝜽𝜽ˆ B (𝑉𝑉) = arg min � 𝐿𝐿M (𝒖𝒖, 𝜽𝜽) 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) 𝑑𝑑𝜽𝜽, 𝒖𝒖

ℬ

(39)

(40)

i.e., the matched average of the emitter positions under the likelihood restricted to ℬ. This estimator is the matched mean of the normalized likelihood whereas the global maximum likelihood estimator of Eq. (30) is its mode, i.e., the two are different functionals of the same normalized likelihood. A network trained on frames drawn from the uniform prior therefore approaches the matched mean of the normalized likelihood whereas EM-GML reaches its mode, so the two differ as the mean and the mode of the same normalized likelihood, which is the origin of the difference between them observed in Section 5. The network approximation of 𝜽𝜽ˆ B is valid only on the prior it was trained on. By the second consequence stated in Section 2.2, the optimality of 𝜽𝜽ˆ B and hence the accuracy of the network hold for frames drawn from 𝑝𝑝(𝜽𝜽), so all comparisons in this paper are made on such frames. 3.2 Network architectures

Five networks are compared. They are not proposed as new architectures for localization and no claim is made that any of them is specific to or optimal for the problem. They are chosen to span a range of inductive biases, so that a result common to all of them can be attributed to learning from frames rather than to the design of any one network. They form four levels. Level 1, learnable activations on a plain network. The adaptive piecewise linear network, denoted APL, applies a linear map followed by a learnable piecewise linear activation built on a grid of knots [31, 32]. It is the oldest and simplest of the five and is included as a deliberate plain baseline, that is, as the network least likely to possess any advantage for this problem. Level 2, learnable edge functions. The Kolmogorov-Arnold network, denoted KAN, replaces the linear map by a learnable univariate function on every edge, each function being the sum of a sigmoid linear unit residual and a B-spline on a fixed grid [33]. KAN is a popularly applied novel network architecture. It is a different and more expressive functional basis than APL at a comparable depth. Level 3, convolutional front end. The networks CNN-APL and CNN-KAN, where CNN stands for convolutional neural network [34], prepend a two layer convolutional stem to the APL and KAN heads respectively. The stem supplies spatial translation aware features that neither plain network can form directly, thus improving training speed and localization accuracy. Level 4, self attention. The vision transformer, denoted ViT [35], treats each pixel of the frame as one token, prepends a learnable class token and applies a stack of transformer encoder blocks with multi-head self attention. The class token output is mapped to the 𝑀𝑀 estimated positions. Self attention captures spatial correlations of a frame in a manner analogous to convolution, whereas it uses no fixed kernel and no built-in translation equivariance. All five networks accept the frame as a vector of 𝐾𝐾𝑥𝑥 𝐾𝐾𝑦𝑦 pixel values normalized by the largest value of the frame and produce 2𝑀𝑀 numbers that are read as 𝑀𝑀 positions in nanometers. The normalization removes the absolute photon scale from the input, so the networks are given the shape of the frame rather than its photon count. The reader is referred to the references for the detailed architectures of the five networks.

3.3 Training All five networks are trained by minimizing Eq. (38) with the matched loss 𝐿𝐿M of Eq. (27) evaluated by the Hungarian assignment. The assignment is recomputed for every frame at every step, so no ordering of the estimates is imposed and the network is free to produce the 𝑀𝑀 positions in any order. Training frames are generated on the fly, that is, a fresh batch of emitter configurations is drawn from the prior and a fresh batch of frames is generated at every step, so no training frame is presented twice and 𝑁𝑁t in Eq. (38) equals the total number of frames generated. The optimizer is Adam [36] with a cosine annealing learning rate [37]. All networks are trained with the same batch size, the same number of steps and therefore the same number of training frames, so that their accuracies and their convergence can be compared on an equal budget. A separate network is trained for each emitter number 𝑀𝑀. 4.

Performance of Bayes estimators for two emitters

The matched Bayes estimator 𝜽𝜽ˆ B of Eq. (28) has no closed form for 𝑀𝑀 ≥ 2 and is the minimizer of the reported metric, so it must be computed before any estimator can be measured against it. We consider the case of 𝑀𝑀 = 2 where the complexity of numerical computations is acceptable. For 𝑀𝑀 = 2 the space of emitter positions 𝜽𝜽 is four dimensional and the posterior can be integrated numerically by quadrature, which yields 𝜽𝜽ˆ B and the degenerate MMSE estimator 𝜽𝜽ˆ MMSE of Eq. (17) exactly up to the grid in the numerical evaluation. This section states the computation, whereas the validation that no estimator falls below 𝜽𝜽ˆ B on frames drawn from the prior is reported in Section 5. 4.1 The posterior on a grid

The prior 𝑝𝑝(𝜽𝜽) is uniform on the central 3 × 3 pixel block ℬ of FOV, both emitters are drawn independently and uniformly from that region, so by Eq. (12) the posterior is proportional to the likelihood and no prior factor enters the computation. The region ℬ is covered by a square grid of 𝑆𝑆 = 𝐺𝐺 2 single emitter nodes 𝝃𝝃1 , … , 𝝃𝝃𝑆𝑆 of spacing Δ𝑔𝑔 where 𝐺𝐺 = 3Δ𝑥𝑥 /Δ𝑔𝑔 with Δ𝑥𝑥 = Δ𝑦𝑦 is the number of nodes in one axis. The posterior being discretized is the joint over the two labeled emitter positions (𝜽𝜽1 , 𝜽𝜽2 ), whose support is the product ℬ × ℬ, i.e. emitter 1 has its own coordinate and emitter 2 has its own. The two emitters are placed at ordered pairs of nodes. A grid point of that product space assigns emitter 1 to node 𝝃𝝃𝑎𝑎 and emitter 2 to node 𝝃𝝃𝑏𝑏 , so (𝝃𝝃𝑎𝑎 , 𝝃𝝃𝑏𝑏 ) and (𝝃𝝃𝑏𝑏 , 𝝃𝝃𝑎𝑎 ) are two distinct points with two different label assignments. There are exactly 𝑆𝑆 2 such ordered pairs. For example, for a 3 × 3 pixel block ℬ, if Δ𝑥𝑥 = Δ𝑦𝑦 = 100 nm, Δ𝑔𝑔 = 6 nm, then 𝐺𝐺 = 50 and there are 𝑆𝑆 = 2,500 nodes for one emitter and 𝑆𝑆 2 = 6,250,000 ordered pairs of nodes for two emitters. For an emitter at node 𝝃𝝃, the signal mean in the 𝒌𝒌th pixel is 𝑠𝑠(𝒌𝒌; 𝝃𝝃) of Eq. (2), so by Eqs. (3) and (6) the mean photon count of the 𝒌𝒌th pixel by the pair (𝝃𝝃𝑎𝑎 , 𝝃𝝃𝑏𝑏 ) is 𝑣𝑣(𝒌𝒌; 𝝃𝝃𝑎𝑎 , 𝝃𝝃𝑏𝑏 ) = 𝑠𝑠(𝒌𝒌; 𝝃𝝃𝑎𝑎 ) + 𝑠𝑠(𝒌𝒌; 𝝃𝝃𝑏𝑏 ) + 𝑏𝑏(𝒌𝒌).

(41)

𝐿𝐿𝑎𝑎,𝑏𝑏 = �[𝑉𝑉(𝒌𝒌) ln 𝑣𝑣 (𝒌𝒌; 𝝃𝝃𝑎𝑎 , 𝝃𝝃𝑏𝑏 ) − 𝑣𝑣(𝒌𝒌; 𝝃𝝃𝑎𝑎 , 𝝃𝝃𝑏𝑏 )] .

(42)

The log posterior of the pair over the grid equals the log-likelihood of Eq. (8) up to an additive constant that is independent of the pair,

𝒌𝒌∈Ω

The normalized posterior weight of the pair is 𝑃𝑃𝑎𝑎,𝑏𝑏 =

𝑒𝑒 𝐿𝐿𝑎𝑎,𝑏𝑏 , ∑𝑆𝑆𝑎𝑎′=1 ∑𝑆𝑆𝑏𝑏′=1 𝑒𝑒 𝐿𝐿𝑎𝑎′,𝑏𝑏′

(43)

which sums to one over the 𝑆𝑆 2 ordered pairs and is exchangeable, i.e., 𝑃𝑃𝑎𝑎,𝑏𝑏 = 𝑃𝑃𝑏𝑏,𝑎𝑎 by the symmetry of Eq. (41) in the two nodes. 4.2 The labeled MMSE estimator

The MMSE estimator of Eq. (17) is the posterior mean of each emitter index and is labeled because it is formed index by index, i.e., the first estimate is the posterior mean of the first emitter position and the second estimate is the posterior mean of the second, before any matching removes the labels. The marginal weight of the first index is 𝑆𝑆 (1) 𝑃𝑃𝑎𝑎 = � 𝑃𝑃𝑎𝑎,𝑏𝑏

(44)

𝑏𝑏=1

and the first estimate is the corresponding weighted mean node 𝑆𝑆

(1) 𝜽𝜽ˆ MMSE,1 = � 𝑃𝑃𝑎𝑎 𝝃𝝃𝑎𝑎 .

(45)

𝑎𝑎=1

The second marginal and the second estimate are formed in the same way from the second index. Since 𝑃𝑃𝑎𝑎,𝑏𝑏 = 𝑃𝑃𝑏𝑏,𝑎𝑎 , the two marginals are identical so the two estimates coincide, which is the degeneracy of Proposition 1. The computed separation �𝜽𝜽ˆ MMSE,1 − 𝜽𝜽ˆ MMSE,2 � is at the level of the floating point precision on every frame, so the degeneracy holds numerically and not only in the mean. 4.3 The matched Bayes estimator The matched Bayes estimator of Eq. (28) minimizes the posterior mean of the matched loss 𝐿𝐿M of Eq. (27). On the grid the objective of a candidate pair (𝒚𝒚1 , 𝒚𝒚2 ) is 𝑆𝑆

𝑆𝑆

𝐽𝐽(𝒚𝒚1 , 𝒚𝒚2 ) = � � 𝑃𝑃𝑎𝑎,𝑏𝑏 min �‖𝒚𝒚𝑖𝑖 − 𝝃𝝃𝑎𝑎 ‖2 + �𝒚𝒚𝑗𝑗 − 𝝃𝝃𝑏𝑏 � , 𝑖𝑖, 𝑗𝑗 = 1,2, 𝑖𝑖 ≠ 𝑗𝑗�

(46)

𝜽𝜽ˆ B = arg min 𝐽𝐽(𝒚𝒚1 , 𝒚𝒚2 ) .

(47)

𝐼𝐼𝑎𝑎,𝑏𝑏 = 𝟏𝟏{‖𝒚𝒚1 − 𝝃𝝃𝑎𝑎 ‖2 + ‖𝒚𝒚2 − 𝝃𝝃𝑏𝑏 ‖2 ≤ ‖𝒚𝒚1 − 𝝃𝝃𝑏𝑏 ‖2 + ‖𝒚𝒚2 − 𝝃𝝃𝑎𝑎 ‖2 }

(48)

𝑎𝑎=1 𝑏𝑏=1

and the estimator is its minimizer

𝒚𝒚1 ,𝒚𝒚2

The minimizer is found by a Lloyd iteration that alternates an assignment step and an update step. The assignment step matches every pair to the candidate by the smaller of the two terms of Eq. (46), which sets by the indicator function 𝟏𝟏(∙) where 𝐼𝐼𝑎𝑎,𝑏𝑏 = 1 selects the identity match and 𝐼𝐼𝑎𝑎,𝑏𝑏 = 0 selects the swap. The update step sets each candidate to the posterior weighted mean of the nodes assigned to it, 𝑆𝑆

𝑆𝑆

𝒚𝒚1 ← � � 𝑃𝑃𝑎𝑎,𝑏𝑏 �𝐼𝐼𝑎𝑎,𝑏𝑏 𝝃𝝃𝑎𝑎 + �1 − 𝐼𝐼𝑎𝑎,𝑏𝑏 � 𝝃𝝃𝑏𝑏 �

(49)

𝒚𝒚2 ← � � 𝑃𝑃𝑎𝑎,𝑏𝑏 �𝐼𝐼𝑎𝑎,𝑏𝑏 𝝃𝝃𝑏𝑏 + �1 − 𝐼𝐼𝑎𝑎,𝑏𝑏 � 𝝃𝝃𝑎𝑎 �

(50)

𝑎𝑎=1 𝑏𝑏=1 𝑆𝑆

𝑆𝑆

𝑎𝑎=1 𝑏𝑏=1

where the weight of each candidate sums to one so no normalization is needed. The two steps do not increase 𝐽𝐽 at any iteration and converge to a fixed point. To keep the fixed point from being a local minimizer the iteration is started from the pair of largest posterior weight, i.e., the

maximum a posteriori pair arg max(𝑎𝑎,𝑏𝑏) 𝑃𝑃𝑎𝑎,𝑏𝑏 , and from several random pairs, and the fixed point of smallest 𝐽𝐽 is kept as 𝜽𝜽ˆ B . By Eqs. (27) and (36) the objective 𝐽𝐽 equals the posterior mean of 𝑀𝑀𝐷𝐷H2 for 𝑀𝑀 = 2 so 𝜽𝜽ˆ B minimizes the reported metric by construction. The iteration uses the posterior 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉)of Eq. (12) and hence the true system parameters, so 𝜽𝜽ˆ B computed this way uses the full system model at every frame, whereas the network of Section 3 approximates the same function of the frame in a single forward pass. 4.4 Grid resolution and Monte Carlo error

The computation is exact up to two controlled errors. The first is the grid spacing Δ𝑔𝑔 whose quantization standard deviation is, by the uniform distribution of an emitter position in the square of grid Δ𝑔𝑔 , 𝜎𝜎𝑔𝑔 =

Δ𝑔𝑔

√12

(51)

A spacing of Δ𝑔𝑔 = 5 to 6 nm gives 𝜎𝜎𝑔𝑔 below 2 nm which is far smaller than the errors of the estimators compared in Section 5 so the grid does not affect the comparison. The second is the Monte Carlo average over the frames of one emitter separation. With 𝑁𝑁 = 200 frames per separation, the relative standard error of the reported ARMSMD-H is near 5 percent and it falls as 𝑁𝑁 −1/2 so the full random dataset of 𝑁𝑁 = 1000 frames reduces it to near 2 percent. 5.

Simulation results

5.1 Simulation configuration System parameters. The camera has 7 × 7 pixels of size Δ𝑥𝑥 = Δ𝑦𝑦 = 100 nm, so the FOV is 700 × 700 nm2, and the frame time is 0.01 s. The PSF is Gaussian with standard deviation 108.81 nm. Every emitter emits at 𝐼𝐼 = 3 × 105 photons per second, i.e., about 3000 photons per emitter per frame. The mixed noise is a frozen non-uniform field, i.e., the background density 𝑏𝑏𝑘𝑘 is a smooth spatial cloud on [4, 6] photons/nm2/s with a correlation length of 2 pixels while the readout density 𝐺𝐺𝑘𝑘 = 𝜇𝜇𝑘𝑘 is drawn independently per pixel and uniformly on [2, 4] with its mean set equal to its variance so the Poisson approximation of Section 2.1 applies. The field is generated once and saved, so every frame and every estimator consumes the identical 𝑏𝑏𝑘𝑘 and 𝐺𝐺𝑘𝑘 that produced the data. The prior draws the 𝑀𝑀 emitter positions independently and uniformly from the central 3 × 3 pixel block ℬ, which is the region on which the networks are trained. Testing datasets. Two datasets are used in testing the trained networks and benchmarking by the EM-GML and UGIA-F estimators. First, Dataset-Random draws a fresh emitter configuration from the prior for each of 𝑁𝑁 = 1000 frames per emitter number 𝑀𝑀, so it is drawn from the same distribution the networks are trained on and it is the in-distribution test. Second, Dataset-Fixed places the 𝑀𝑀 emitters equally spaced on a circle of radius 𝑟𝑟, so the nearestneighbor separation is 2𝑟𝑟 sin(𝜋𝜋/𝑀𝑀), which decreases as 𝑀𝑀 grows at fixed 𝑟𝑟. Each pair (𝑟𝑟, 𝑀𝑀) is repeated over 𝑁𝑁 = 25 noise realizations with 𝑟𝑟 ∈ {25, 50, 75, 100, 125, 150} nm and 𝑀𝑀 ∈ {1, … , 5}. The circle constellations are not drawn from the training prior, so Dataset-Fixed probes the generalization of the networks to configurations they were not trained on. Network parameters. The five networks are sized to a common scale so that a result common to all of them reflects the inductive bias rather than the capacity. APL applies a linear map to two hidden layers of 256 and 128 units where every unit carries a learnable piecewise linear activation on 16 knots spaced uniformly on [−3, 3]. KAN uses two hidden layers of 64 and 32 units where every edge is a learnable univariate function that sums a sigmoid linear unit residual and a cubic B-spline on a grid of 5 intervals, i.e., a spline order of 3. CNN-APL and CNN-KAN prepend a two-layer convolutional stem of 16 channels with a 3 × 3 kernel and unit padding to the APL and KAN heads respectively, so the stem maps the 7 × 7 frame to 16 15

feature maps that flatten to 7 × 7 × 16 = 784 values feeding the head, where the APL head carries two hidden layers of 256 units and the KAN head carries hidden layers of 64 and 32 units. ViT reads each of the 49 pixels as a token of a 128 dimensional embedding, prepends a learnable class token and applies 6 transformer encoder layers of 4 self-attention heads each with a feedforward width of 512 and a dropout of 0.1, then maps the class token through a twolayer head of width 512 to the output, which totals about 1.25 million weights. Every network ends in a layer of 2𝑀𝑀 outputs read as 𝑀𝑀 positions in nanometers. Training. All five networks are trained by minimizing the matched Hungarian loss 𝐿𝐿M of Eq. (27) on frames generated on the fly from the prior, with the Adam optimizer and a cosine annealing learning rate over the run. Every network is trained with a batch size of 256 for 40000 steps, so each network sees 1.024 × 107 training frames, and a separate network is trained for each 𝑀𝑀. The learning rate is annealed from 10−3 to 10−5 for APL, CNN-APL, KAN and CNN-KAN and from 10−4 to 10−6 for ViT. Evaluation. Every estimator is evaluated by the ARMSMD-H of Eq. (37), i.e., the root mean square minimum distance under the Hungarian assignment averaged over the frames. The UGIA-F and EM-GML estimators use the true system parameters at test time and, for UGIAF, the true positions in addition, whereas the networks use the system parameters only in training to synthesize the frames and receive only the frame at test time. The EM-GML runs with multiple random initial positions in the support ℬ to reach the global maximizer reliably, whereas a fully practical run initializes the EM by a rough estimator such as SIC [21] that uses no system information, so EM-GML is a practical estimator. Python code. We developed the SMLM_Lib, a Python library for SMLM [38]. Based on it, custom Python code was developed for the simulations in this paper. The Python code can detect and take use of a GPU if it is available on a computer. All simulations were carried out on a Dell workstation XPS-8960 with Intel Core i7 of 2.10 GHz, 64 GB RAM and NVIDIA GeForce RTX 4070 of 12 GB VRAM. The training times are listed in Table 1. The CNN head can speed up the convergence of APL and KAN in the training. Six example training frames are shown in Fig. 1. Table 1. Training time (min:sec) with 10.24 million training frames. CNN-

KAN

CNN-

ViT

𝑀𝑀

APL

4:43

2:19

7:16

3:58

13.14

5:11

2:22

7:18

4:03

13:05

4:50

2:19

7:17

4:00

13:07

4:42

2:33

7:23

3:59

13:07

4:42

2:21

7:17

3:56

13:06

APL

KAN

Fig. 1. Six example training frames. Top row 𝑀𝑀 = 0, 1, 2 and bottom row 𝑀𝑀 = 3, 4, 5, where 𝑀𝑀 = 0 is a pure-noise frame showing the frozen non-uniform noise field. True emitter positions are shown as red dots. With the area of 0.09 µm2 for the central 3 × 3 pixel block ℬ, the emitter density is 22.2, 33.3, 44.4, 55.6 emitters/µm2 for 𝑀𝑀 = 2, 3, 4, 5 respectively.

On the same workstation a trained network estimates one frame in about 0.5 to 16 µs at a saturating batch, i.e., between 6 × 104 and 2 × 106 frames per second or across the five networks. The multi-start EM of EM-GML at 𝐾𝐾 = 250 starts and 300 iterations takes about 0.5 seconds per frame, so the networks are four to six orders of magnitude faster than the multistart EM. 5.2 Two emitters against the Bayes optimum

For 𝑀𝑀 = 2 the optimum matched Bayes estimator 𝜽𝜽ˆ B and the MMSE estimator 𝜽𝜽ˆ MMSE are computed by the quadrature of Section 4, so the estimators can be measured against the Bayes optimum on the same frames. Table 2 reports the accuracy and the agreement of every estimator on the 1000 frames of Dataset-Random, i.e., on frames drawn from the prior. The accuracy is the ARMSMD-H against the true positions, whereas the agreement is the ARMSMD-H against the Bayes estimate on the same frame, i.e., the direct measure of how closely an estimator reproduces 𝜽𝜽ˆ B frame by frame. Estimator Accuracy Agreement

Table 2. Accuracy and Agreement with 𝜽𝜽ˆ 𝐁𝐁 in nm for 𝑴𝑴 = 𝟐𝟐.

Bayes

EM-

𝜽𝜽ˆ B 13.21

GML

ViT

CNN-

KAN

APL

KAN

APL

UGIAF

MMSE

14.03

14.37

14.53

14.90

16.60

17.18

17.19

85.97

4.64

5.25

6.25

7.46

10.31

11.05

21.53

84.96

Three results follow from Table 2. First, no estimator falls below the Bayes optimum, i.e., 𝜽𝜽ˆ B attains the smallest ARMSMD-H of 13.21 nm and a paired test over the frames confirms that no estimator is significantly below it, which validates the whole pipeline against the bound of Section 2.5. Second, every network is at or below UGIA-F, and the three stronger networks ViT, CNN-KAN and CNN-APL beat UGIA-F and approach both EM-GML and the Bayes optimum, with ViT within 1.2 nm of 𝜽𝜽ˆ B and essentially equal to EM-GML. This holds although the networks receive only the frame at test time. Third, the MMSE estimator is degenerate as Proposition 1 predicts, i.e., its two estimates coincide to 7 × 10−5 nm on every 17

frame, which is the residual of the finite-precision quadrature rather than a nonzero separation, and its ARMSMD-H of 85.97 nm is that of the centroid. The agreement row separates the networks from UGIA-F. EM-GML and ViT agree with the Bayes estimate to 4.64 and 5.25 nm, i.e., they reproduce 𝜽𝜽ˆ B frame by frame and not only on average, whereas UGIA-F agrees only to 21.53 nm although its accuracy is similar to that of APL, i.e., it reaches a comparable average error by a structurally different route. The distinction is confirmed by the excess risk, i.e., the squared distance of an estimator from the true positions decomposes almost exactly into its squared distance from the Bayes estimate plus the squared distance of the Bayes estimate from the truth for the networks and for EM-GML. For example, for ViT, 14.372 ≈ 5.252 + 13.212 . In contrast, the decomposition fails for UGIA-F. This supports reading the small excess of a network over the Bayes optimum as approximation error rather than a different kind of error. The dependence on the emitter separation is shown in Fig. 2, i.e., the ARMSMD-H of every estimator against the separation 𝑑𝑑 at fixed 𝑑𝑑. UGIA-F diverges as the two emitters merge, i.e., it rises from about 10 nm at 𝑑𝑑 = 250 nm to 116 nm at 𝑑𝑑 = 10 nm as the Fisher information becomes near singular, whereas 𝜽𝜽ˆ B and the networks remain bounded near 18 nm. The curves cross, i.e., UGIA-F is competitive only at large separation while the networks and 𝜽𝜽ˆ B dominate at small separation. Fig. 2 conditions on a fixed separation, so it is a conditional view of the prior, and an estimator may fall below 𝜽𝜽ˆ B at a single separation without contradicting the optimality of Section 2.5, which is an average over the prior. The optimality is the prioraveraged statement of Table 2, whereas Fig. 2 shows where in the configuration space each estimator gains or loses.

Fig. 2. Accuracy of estimators versus the Bayes optimum against the separation 𝑑𝑑 for 𝑀𝑀 = 2 emitters.

5.3 Accuracy versus emitter number on the Random dataset

Table 3 reports the ARMSMD-H of every estimator on Dataset-Random for 𝑀𝑀 = 1 to 5, i.e., the in-distribution accuracy as the number of emitters grows. With the area of 0.09 µm2 for the central 3 × 3 pixel block ℬ where the emitters are distributed, the emitter density is 22.2, 33.3, 44.4, 55.6 emitters/µm2 for 𝑀𝑀 = 2, 3, 4, 5 respectively, which are quite high. All estimators are evaluated on the same 1000 frames per 𝑀𝑀 and the UGIA-F Gaussian draw is seeded so the run is reproducible. Fig. 3 shows example images of estimated positions by the seven estimators for 𝑀𝑀 = 5. More example images for 𝑀𝑀 = 1 to 4 can be found in Ref. [39].

Table 3. ARMSMD-H (nm) on Dataset-Random. CNN-

KAN

CNN-

ViT

UGIA-

EM-

GML

𝑀𝑀

APL

9.63

9.34

9.35

9.77

9.46

9.49

17.22

14.90

16.60

14.53

14.37

25.80

14.02

29.53

20.71

24.63

22.56

21.27

49.65

20.61

36.30

28.02

32.11

28.55

27.26

116.35

26.11

37.73

31.67

35.64

33.22

31.59

172.39

30.30

APL

KAN

Fig. 3. Example images of positions estimated by the seven estimators for 𝑀𝑀 = 5. The RMSMD-H for each image is shown.

For a single emitter there is no overlap and every estimator attains the same ARMSMD-H near 9.5 nm. As 𝑀𝑀 grows, the emitters overlap more often and UGIA-F degrades steeply, i.e., it rises to 172 nm at 𝑀𝑀 = 5 because the CRB diverges whenever a pair of emitters merges. EM-GML remains the most accurate over the whole range and stays bounded, i.e., it grows only from 9.5 to 30 nm. Every network stays far below UGIA-F for 𝑀𝑀 ≥ 3 and approaches EM-GML, i.e., at 𝑀𝑀 = 5 the best network ViT attains 31.6 nm against 30.3 nm for EM-GML and 172 nm for UGIA-F. The networks order themselves by architecture, i.e., the convolutional and attention networks ViT, CNN-KAN and CNN-APL lead while the plain APL and KAN trail, which shows that the advantage comes from learning from frames rather than from any single design. The hypothesis of the paper is therefore confirmed on the in-distribution test, i.e., the networks outperform the UGIA-F benchmark and approach EM-GML. UGIA-F is realized once per frame from Eq. (34), i.e., a single Monte Carlo draw on the same footing as the single noisy frame that every other estimator receives, so it is reported from one realization rather than averaged. The aggregate is a high-variance quantity at overlap, i.e., it is dominated by the rare near-coincident frames where the CRB diverges. As sown in Table 5 of Section 5.5, for 𝑀𝑀 = 2 the 17 frames with 𝑑𝑑min < 25 nm carry a conditional UGIA-F error of 162 nm, so a handful of frames lift the movie-wide value. The aggregate therefore varies between realizations, e.g., 25.8 nm here and 17.2 nm in the validation run of Table 2 on the same configurations, which is the sampling spread of the oracle rather than a discrepancy. 5.4 Generalization to the circle constellations

Dataset-Fixed places the emitters on a circle of radius 𝑟𝑟, so the constellations are not drawn from the training prior and the nearest-neighbor separation 2𝑟𝑟 sin(𝜋𝜋/𝑀𝑀) tightens as 𝑀𝑀 grows. This dataset therefore tests whether the networks generalize beyond the configurations they were trained on. The number of Monte Carlo runs is 𝑁𝑁 = 25. Table 4 reports the ARMSMDH of every estimator for each radius 𝑟𝑟 and emitter number 𝑀𝑀. Fig. 4 shows the estimated 𝑁𝑁 = 25 positions by the seven estimators in comparison of the true positions for 𝑟𝑟 = 125 nm and 𝑀𝑀 = 5. All images of estimated positions by the seven estimators for 𝑟𝑟 = 25, 50, 75, 100, 125 nm and 𝑀𝑀 = 1, 2, 3, 4, 5 can be found in Ref. [39]. Table 4. ARMSMD-H (nm) on Dataset-Fixed.

𝑟𝑟

CNN-

UGIA-

EM-

GML

9.49

9.18

18.53

21.96

22.14

19.31

159.14

22.65

15.51

15.06

1625.92

23.77

19.93

13.62

280.29

21.34

6.85

6.93

10.79

6.77

15.61

15.43

13.50

15.54

30.10

30.74

49.12

30.93

27.42

28.18

27.40

185.83

31.78

20.87

29.01

22.34

1205.47

34.73

APL

𝑀𝑀

9.39

9.09

9.27

9.11

9.46

19.02

19.55

18.19

19.86

22.42

21.15

21.21

18.71

20.76

19.14

16.82

23.55

19.86

15.44

6.85

6.58

6.74

17.89

17.28

15.90

33.38

33.95

30.88

26.78

29.14

21.53

26.21

8.30

8.32

8.16

8.39

7.98

9.46

8.29

12.32

10.96

11.88

11.10

11.35

11.20

10.41

31.75

34.54

31.81

27.69

27.87

23.64

25.22

29.32

36.07

29.83

36.82

37.38

53.58

36.61

31.30

37.72

27.69

38.67

40.46

252.28

41.91

100

9.93

9.55

9.39

9.25

10.01

9.03

9.40

100

10.45

9.60

9.58

8.92

10.43

10.64

9.27

100

15.74

20.69

17.70

14.83

16.78

12.82

13.71

100

34.76

30.69

29.33

33.67

41.14

39.32

32.59

100

38.81

48.73

32.05

46.70

47.39

71.73

45.35

125

8.92

9.12

9.07

9.14

9.42

11.58

9.10

125

10.41

10.30

10.49

10.48

10.96

10.98

10.03

125

14.06

12.02

13.86

12.59

12.50

12.10

11.80

125

29.09

22.19

23.25

21.09

30.02

22.82

24.33

125

42.56

57.36

28.95

56.07

50.87

50.14

34.62

150

9.79

9.00

9.22

9.11

9.55

8.83

9.79

150

8.33

9.20

9.42

8.82

9.52

9.58

8.79

150

12.22

12.74

18.36

12.82

12.89

10.16

10.80

150

29.71

15.46

17.41

20.16

18.21

16.21

15.90

150

53.51

74.56

31.64

66.98

57.15

27.20

35.32

(nm)

APL

KAN

ViT

Fig. 4. The images of positions estimated by the seven estimators over 𝑁𝑁 = 25 frames for a radius of 𝑟𝑟 = 125 nm and 𝑀𝑀 = 5 emitters.

The dominant feature of Table 4 is the divergence of UGIA-F at severe overlap, i.e., the smallest radius with the largest emitter number packs the circle within the PSF so the Fisher information becomes near singular and the best unbiased accuracy is unbounded, e.g., UGIAF reaches 1626 nm at 𝑟𝑟 = 25 nm and 𝑀𝑀 = 4 , and 1205 nm at 𝑟𝑟 = 50 nm and 𝑀𝑀 = 5 , whereas no learned or likelihood estimator exceeds 75 nm anywhere in the table. This is the CRB divergence of Section 2.7 made concrete, i.e., the biased EM-GML and the networks remain bounded where the oracle UGIA-F does not. At large radius all estimators converge to the single-emitter ARMSMD-H near 9 to 15 nm, i.e., the overlap is resolved and the estimators agree. Between these limits the networks generalize well but not perfectly. They track EM-GML across most of the table, whereas at the tightest packings and the largest 𝑀𝑀 some networks exceed EM-GML, e.g., CNN-APL reaches 74.6 nm at 𝑟𝑟 = 150 nm and 𝑀𝑀 = 5 . This is consistent with Section 3.1, i.e., the network approximation of the Bayes estimator holds on the training prior and degrades on constellations far from it, so Dataset-Fixed marks the edge of generalization while Dataset-Random of Section 5.3 is the fair in-distribution comparison. 5.5 Stratified analysis by emitter overlap The aggregate of Section 5.3 averages over all overlap severities, so it hides where in the configuration space a network overtakes a benchmark. The per-frame results of DatasetRandom are therefore stratified by the minimum emitter separation 𝑑𝑑min of each frame, which isolates the effect of overlap from the effect of train-test mismatch because every stratum is drawn from the training prior. Table 5 reports the conditional ARMSMD-H within each stratum for 𝑀𝑀 = 2, 3, 4, 5 . The case 𝑀𝑀 = 1 is excluded because a single emitter has no pairwise separation and 𝑑𝑑min is undefined. Table 5. Conditional ARMSMD-H on Dataset-Random by minimum emitter separation 𝒅𝒅𝐦𝐦𝐦𝐦𝐦𝐦 . 𝑀𝑀

𝑑𝑑𝑚𝑚𝑚𝑚𝑚𝑚

APL

CNN-

KAN

CNN-

ViT

UGIA-

EM-

GML

[0, 25)

𝑛𝑛

22.71

20.19

19.87

20.13

20.78

161.89

26.44

[25, 50)

18.19

15.69

16.24

16.88

16.15

29.07

20.80

[50, 75)

19.11

17.82

19.25

18.59

19.24

24.34

18.53

(nm)

APL

KAN

[75, 100)

114

21.05

17.60

19.82

18.49

17.85

17.80

16.75

[100, 150)

231

18.99

15.35

18.14

14.81

15.43

12.77

13.68

≥150

516

14.68

13.22

14.51

12.11

11.51

10.85

11.03

[0, 25)

29.28

19.83

25.92

24.95

21.28

156.39

23.16

[25, 50)

180

27.33

20.48

23.22

21.60

21.03

47.38

22.50

[50, 75)

182

28.06

20.89

23.97

21.34

21.58

30.31

23.38

[75, 100)

191

28.39

22.29

27.09

25.01

22.63

21.17

21.40

[100, 150)

261

32.73

20.80

26.12

23.72

22.20

14.35

18.54

≥150

115

29.45

18.19

18.42

16.55

15.99

12.39

12.86

[0, 25)

116

33.31

26.93

28.94

24.73

24.07

307.99

25.86

[25, 50)

264

36.25

27.38

32.65

28.18

26.87

85.46

25.71

[50, 75)

263

36.15

27.95

31.92

28.62

26.95

36.39

27.69

[75, 100)

214

36.09

28.37

31.73

29.50

28.54

30.49

26.35

[100, 150)

134

39.50

30.18

34.75

30.83

29.55

20.10

23.95

≥150

34.37

19.26

28.23

24.23

17.61

19.98

16.37

[0, 25)

199

36.95

29.81

34.17

31.76

30.16

354.44

29.05

[25, 50)

349

36.98

31.59

34.67

31.42

29.99

104.74

29.44

[50, 75)

298

37.84

31.15

36.04

33.18

32.62

48.29

29.96

[75, 100)

123

39.33

34.01

38.04

37.34

34.05

37.88

34.20

[100, 150)

43.15

38.76

41.57

43.43

37.50

24.86

34.19

Note: The count 𝑛𝑛 is the number of frames in the stratum. A network value is set in bold where it is smaller than the EM-GML value in the same stratum and the same 𝑀𝑀, whereas the oracle UGIA-F is never marked.

Two patterns stand out. First, the advantage of the networks over the oracle UGIA-F is concentrated in the overlapped strata, i.e., at the smallest 𝑑𝑑min the networks are far more accurate than the diverging UGIA-F, e.g., for 𝑀𝑀 = 2 the best network reaches 19.9 nm against 162 nm for UGIA-F in the [0, 25) stratum, whereas at 𝑑𝑑min ≥ 150 nm the oracle is best and the networks trail by a few nanometers. Every network is more accurate than UGIA-F throughout the overlapped strata below 75 nm at all 𝑀𝑀. Second, the networks also beat EMGML in the most overlapped strata at low emitter number, which the bold entries of Table 5 mark, i.e., for 𝑀𝑀 = 2 every network is more accurate than EM-GML in the two strata below 50 nm, whereas for 𝑀𝑀 = 5 EM-GML leads in every stratum except that ViT and CNN-APL edge past it in the [75, 100) stratum. 5.6 Verification of the EM-GML benchmark

EM-GML serves as the maximum likelihood benchmark in every table, so the multi-start EM that realizes it must be shown to reach the global maximum of the likelihood rather than a local one. For 𝑀𝑀 = 2 the parameter space of emitter positions is four dimensional and can be searched exhaustively, which provides a reference that the multi-start is measured against. An exhaustive grid search over the whole parameter space at a grid step of 5 nm gives the maximizing pair of grid positions on each frame, which is then refined by EM to the exact offgrid maximum of the mode it occupies and is taken as the reference global maximum. On a frame the multi-start is counted as reaching the reference when its largest log-likelihood over the 𝐾𝐾 starts is within a tolerance of the reference log-likelihood. The comparison is on the loglikelihood rather than on the position because two configurations far apart in position can be nearly tied in likelihood when the surface is flat.

Fig. 5 shows the fraction of the 1000 frames on which the multi-start reaches the reference and the mean log-likelihood shortfall below it against the number of starts 𝐾𝐾. The fraction rises from 66 percent at 𝐾𝐾 = 1 to 99.7 percent at 𝐾𝐾 = 250 whereas the mean shortfall falls by more than three orders of magnitude to about 10−5 , i.e., even on the few frames that do not reach the reference exactly the remaining gap is negligible. A separate check finds that raising the EM iteration budget from 300 to 1000 changes the fraction little, so the budget of 300 iterations used throughout is sufficient. The multi-start EM with 𝐾𝐾 = 250 therefore attains the global maximum likelihood estimate reliably for 𝑀𝑀 = 2, which justifies the use of EM-GML as the GML benchmark in this paper.

Fig. 5. Verification of the EM-GML benchmark against an exhaustive grid search for 𝑀𝑀 = 2. The navy curve is the percentage of the 1000 Dataset-Random frames on which multi-start EM reaches the reference global maximum, whereas the crimson curve is the mean log-likelihood shortfall below the reference, both against the number of EM starts 𝐾𝐾.

Discussion

6.1 Information used by the three estimators and their practical usefulness The three estimators differ in the information they use when they estimate the emitter positions from a given frame. A network uses only the frame, since the system parameters and the true positions enter only in training through the synthesized frames and their labels. EM-GML uses the system parameters, i.e., it maximizes the per-frame likelihood and does not use the true positions. UGIA-F uses the system parameters and the true positions, since its Gaussian draw of Eq. (34) is centered at the true positions. Only UGIA-F therefore requires information that a real experiment cannot supply, so UGIA-F alone is an oracle whereas the networks and EMGML are practically useful estimators. At estimation time the networks use the least, i.e., the frame alone, and estimate by a single non-iterative forward pass that is much faster than the iterative EM of EM-GML as reported in Section 5.1, so once trained they are the most readily deployed. EM-GML needs the calibrated system model of every frame but not the position truth and is deployable wherever the model is known. UGIA-F cannot be run on data whose positions are unknown and serves only as the benchmark of the best accuracy attainable by an unbiased estimator. The comparison is nonetheless made on the same system model, i.e., EM-GML and UGIAF use it directly while the networks embed it in their training frames, and the networks additionally use the prior over configurations. A frame-trained network converts the observed frame into an estimate through the fixed map learned in training, which encodes both the system model of the synthesized frames and the prior, whereas EM-GML maximizes the per-frame likelihood under a flat prior. The prior is the side information that lets a biased estimator stay

bounded where the unbiased CRB diverges, so the networks and EM-GML remain accurate at overlap while UGIA-F does not, and the networks are free to match EM-GML and to outperform UGIA-F on the frames drawn from the prior they were trained on. This is why the results of Section 5.3 and Section 5.5 are not paradoxical. The networks beat the oracle UGIA-F over the overlapped strata not by using more of the system but by carrying the prior, which UGIA-F does not use, so a biased learned estimator stays accurate wherever the prior is informative about the overlapped configurations, which is exactly where the unbiased UGIA-F is weakest. 6.2 The network as a practical approximation of the Bayes estimator Section 2.5 and Section 4 establish that the matched Bayes estimator 𝜽𝜽ˆ B has no closed form for 𝑀𝑀 ≥ 2 and is computed by a fixed-point iteration that uses the posterior of Eq. (12) and hence the full system model. A direct computation of 𝜽𝜽ˆ B therefore uses the posterior and hence the full system model at every frame, i.e., the same per-frame model computation that EM-GML performs, though unlike UGIA-F it does not need the true positions. The network amortizes this computation, i.e., it is trained to approximate the same function of the frame that 𝜽𝜽ˆ B computes and, once trained, evaluates it in a single forward pass without the per-frame posterior integration. The forward pass is non-iterative, so once trained the network estimates far faster than the iterative EM of EM-GML and the iterative fixed point that computes 𝜽𝜽ˆ B , which is an important practical advantage of the network. The simulation supports this reading in two ways. On Dataset-Random of Section 5.3 the networks track EM-GML and beat the oracle UGIA-F over the overlapped strata, which is the behavior expected of an estimator that minimizes the mean square error over the prior. On Dataset-Fixed of Section 5.4 the accuracy degrades on the tightest circle constellations, i.e., exactly the configurations that lie farthest from the prior, which confirms that the network approximates 𝜽𝜽ˆ B on the training prior rather than a parameter-free estimator that would hold everywhere. The scope of the approximation is therefore the prior, and within that scope the network reaches the accuracy of EM-GML, which is the hypothesis of the paper. Hence, the results present a strong promise of neural networks in developing highthroughput large-FOV super spatiotemporal resolution SMLM [40]. 6.3 Necessity of the Hungarian assignment Sections 2.3 to 2.5 theoretically analyze the performance of the MMSE estimator, the sorted Bayes estimator, and the matched Bayes estimator. Correspondingly, their loss functions require no match, one coordinate match, and full match between the estimated position and the true positions. In simulations, Section 5.2 reports the accuracy of the MMSE estimator that generates a single estimated position for both of the 𝑀𝑀 = 2 emitters, confirming the theoretical result in Section 2.3 that the estimated position coordinates in both 𝑥𝑥 and 𝑦𝑦 axes merge to a single position. In addition, the networks in all reported simulation results are trained by using the Hungarian assignment aiming to approach the matched Bayes estimator. Their reconstructed images of estimated emitter positions show balanced small spread over both 𝑥𝑥 and 𝑦𝑦 coordinates, confirming the theoretical analysis in Section 2.5. Though not reported in detail, we also trained the five networks to approach the sorted Bayes estimator by matching only 𝑥𝑥 coordinates between the estimated and the true positions. The result clearly shows that when the 𝑥𝑥 coordinates of two emitter positions are near, their network estimated 𝑦𝑦 coordinates are widely spread, confirming the theoretical analysis in Section 2.4. Hence, for a neural network to approach the optimum matched Bayes estimator, the Hungarian assignment must be applied in training. 6.4 Why the oracle UGIA-F diverges while the biased estimators stay bounded

The most visible feature of Table 3 and Table 4 is the divergence of UGIA-F at severe overlap against the bounded error of EM-GML and the networks. This follows from the status of unbiasedness as a constraint rather than an optimality property, which is stated in Section 2.7. When a pair of emitters merges, the Fisher information matrix becomes near singular and the CRB diverges, so the best accuracy attainable under the unbiasedness constraint is unbounded and UGIA-F inherits that divergence through Eq. (34). EM-GML carries a bias and is therefore free of the constraint, so its error stays bounded where the CRB does not, e.g., it grows only from 9.5 to 30 nm across 𝑀𝑀 = 1 to 5 on Dataset-Random whereas UGIA-F reaches 172 nm. The networks share the bounded behavior of EM-GML rather than the divergence of UGIAF. A network minimizes the matched loss of Eq. (27) averaged over the prior, i.e., it minimizes a mean square error and not a variance under an unbiasedness constraint, so it too is a biased estimator and is not bounded by the CRB. The prior supplies the information that a biased estimator needs to stay accurate when the frame alone is nearly uninformative, i.e., when two emitters overlap the prior still favors the configurations it was trained on and the network settles on a plausible pair rather than on the diverging unbiased solution. The divergence of UGIA-F and the boundedness of the learned estimators are therefore two sides of the same fact, i.e., an unbiased estimator pays for the merging of the emitters through the CRB whereas a biased estimator does not. 7.

Conclusion

This paper tested the hypothesis that a neural network trained on synthesized frames can outperform the unbiased oracle UGIA-F and approach the likelihood estimator EM-GML in multi-emitter localization. Five networks were studied, i.e., APL, KAN, CNN-APL, CNNKAN and ViT, each trained on frames drawn on the fly from the prior and each receiving only the frame at test time, whereas UGIA-F and EM-GML receive the PSF, the emitter intensity and the per-pixel noise densities, and UGIA-F receives the true positions as well. The matched Bayes estimator was identified as the target that a frame-trained network approximates, and Section 4 verified by quadrature on two emitters that this estimator has no closed form and is computed by a fixed-point iteration that needs the full system model, so a direct computation itself uses the full system model at every frame. The simulation confirmed the hypothesis on the in-distribution test. On Dataset-Random every network stayed far below UGIA-F for 𝑀𝑀 ≥ 3 and approached EM-GML, e.g., at 𝑀𝑀 = 5 the best network ViT attained 31.6 nm against 30.3 nm for EM-GML and 172 nm for UGIAF. The networks ordered themselves by architecture while all of them beat the oracle UGIA-F, which shows that the advantage comes from learning from frames rather than from any single design. The stratified analysis located the advantage in the overlapped configurations, i.e., where the emitters merge, the networks were far more accurate than the diverging UGIA-F and, at low emitter number, more accurate than EM-GML as well. Two findings frame the result. First, the divergence of the oracle UGIA-F at overlap and the bounded error of the networks are two sides of unbiasedness being a constraint rather than an optimality property, i.e., a biased estimator is free of the CRB that the merging emitters make unbounded. Second, at estimation time the networks use only the frame, EM-GML uses the system parameters, and UGIA-F uses the system parameters and the true positions, so only UGIA-F is an oracle while the networks and EM-GML are practically useful, and the networks additionally use the prior in training, so their accuracy comes from learning the Bayes-optimal estimator rather than from any information advantage. The approximation holds on the prior it was trained on. On Dataset-Fixed the accuracy degraded on the tightest circle constellations, i.e., the configurations that lie farthest from the training prior, which marks the edge of generalization and confirms that the network approximates the matched Bayes estimator on the prior rather than a parameter-free estimator that would hold everywhere. Within that scope the networks reached the accuracy of EM-GML, so a frame-trained network is a practical estimator for multi-emitter localization, i.e., once

trained on frames synthesized for the imaging conditions of interest it evaluates the Bayesoptimal estimate in a single non-iterative forward pass, which is much faster than the iterative EM that EM-GML requires, presenting a strong promise of neural networks in high-throughput large-FOV SMLM. This paper is a proof of concept. The capability of a network to approach the optimal matched Bayes estimator, i.e., the estimator that minimizes the mean square error over the prior, is established here in concept by simulation on a FOV of 700 × 700 nm2. With the emitters distributed in the central 3 × 3 pixel block of area 0.09 µm2, the corresponding emitter density of 22.2, 33.3, 44.4, 55.6 emitters/µm2 for 𝑀𝑀 = 2, 3, 4, 5 is quite high. This study suggests that supported with an advanced GPU, the convolutional and attention networks CNN-APL, CNNKAN and ViT can be trained on frames of a practically large size, e.g., 2048 × 2048 pixels, to approach the optimum Bayes estimator to eventually realize high-throughput large-FOV super spatiotemporal resolution SMLM. Appendix A: Proofs A.1 Proof of Proposition 1 By Eqs. (2), (3), (5) and (6) the frame mean 𝑣𝑣(𝒌𝒌) depends on 𝜽𝜽 only through the sum 𝑄𝑄(𝒌𝒌), which is symmetric in the emitter index, so 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝐏𝐏𝜋𝜋 𝜽𝜽) = 𝑓𝑓𝑉𝑉 (𝑉𝑉; 𝜽𝜽) for every 𝜋𝜋 ∈ 𝒮𝒮𝑀𝑀 . Together with the exchangeability of the prior and Eq. (12) this gives 𝑝𝑝(𝐏𝐏𝜋𝜋 𝜽𝜽 ∣ 𝑉𝑉) = 𝑝𝑝(𝜽𝜽 ∣ 𝑉𝑉), so the posterior is exchangeable and all its marginals are identical. The conditional means are therefore equal. █

A.2 Proof of Corollary 1

ˆ = 𝒄𝒄 + 𝒆𝒆 and note that 𝒄𝒄 − 𝜽𝜽1 = −(𝒄𝒄 − 𝜽𝜽2 ) and ‖𝒄𝒄 − 𝜽𝜽𝑚𝑚 ‖ = 𝑑𝑑/2 for 𝑚𝑚 = 1,2 . Write 𝒎𝒎 Expanding both squared norms and summing, the cross terms cancel and Eq. (19) follows. █ A.3 Proof of Proposition 2

Take 𝑥𝑥1 ≤ 𝑥𝑥2 without loss of generality and let 𝐴𝐴 = {𝑋𝑋1 < 𝑋𝑋2 } be the event that the order is preserved. The order is reversed exactly when 𝑋𝑋1 > 𝑋𝑋2 , so by Eq. (22) P(𝐴𝐴c ) = 𝑝𝑝 and P(𝐴𝐴) = 1 − 𝑝𝑝. The indicator 𝟏𝟏𝐴𝐴 depends only on 𝜉𝜉1 and 𝜉𝜉2 and is therefore independent of 𝑌𝑌1 and 𝑌𝑌2 . Consider first 𝑌𝑌(1) . Since 𝑌𝑌(1) = 𝟏𝟏𝐴𝐴 𝑌𝑌1 + (1 − 𝟏𝟏𝐴𝐴 )𝑌𝑌2 with 𝟏𝟏𝐴𝐴 independent of (𝑌𝑌1 , 𝑌𝑌2 ), 𝔼𝔼�𝑌𝑌(1) � = (1 − 𝑝𝑝) 𝑦𝑦1 + 𝑝𝑝 𝑦𝑦2 ,

(52)

2 𝔼𝔼�𝑌𝑌(1) � = 𝜎𝜎 2 + (1 − 𝑝𝑝) 𝑦𝑦12 + 𝑝𝑝 𝑦𝑦22 .

(53)

Var�𝑌𝑌(1) � = 𝜎𝜎 2 + 𝑝𝑝(1 − 𝑝𝑝)(𝑦𝑦1 − 𝑦𝑦2 )2 ,

(54)

Subtracting the square of the mean in Eq. (52) from Eq. (53),

which is Eq. (23). Consider next 𝑋𝑋(1) = min(𝑋𝑋1 , 𝑋𝑋2 ) = 𝑆𝑆 − |𝑇𝑇|⁄2 with 𝑆𝑆 = (𝑋𝑋1 + 𝑋𝑋2 )⁄2 and 𝑇𝑇 = 𝑋𝑋2 − 𝑋𝑋1 . Since 𝑋𝑋1 and 𝑋𝑋2 are independent with the common variance 𝜎𝜎 2 , 𝑆𝑆 and 𝑇𝑇 are uncorrelated and jointly Gaussian, hence independent, with Var(𝑆𝑆) = 𝜎𝜎 2 /2 and 𝑇𝑇 ∼ 𝒩𝒩(𝑑𝑑𝑥𝑥 , 2𝜎𝜎 2 ). Therefore 𝜎𝜎 2 1 + [𝔼𝔼(𝑇𝑇 2 ) − 𝔼𝔼2 (|𝑇𝑇|)] 2 4 1 = 𝜎𝜎 2 + [𝑑𝑑𝑥𝑥2 − 𝔼𝔼2 (|𝑇𝑇|)], 4

Var�𝑋𝑋(1) � =

(55)

using 𝔼𝔼[𝑇𝑇 2 ] = 2𝜎𝜎 2 + 𝑑𝑑𝑥𝑥2 . The first absolute moment of the folded normal has the closed form

𝔼𝔼(|𝑇𝑇|) =

2𝜎𝜎

√𝜋𝜋

exp �−

𝑑𝑑𝑥𝑥2 � + 𝑑𝑑𝑥𝑥 (1 − 2𝑝𝑝), 4𝜎𝜎 2

(56)

where 𝜏𝜏 2 = 2𝜎𝜎 2 and 2Φ(𝑑𝑑𝑥𝑥 /𝜏𝜏) − 1 = 1 − 2𝑝𝑝 by Eq. (22). Since 𝑇𝑇 = 𝑋𝑋2 − 𝑋𝑋1 , Eq. (56) is Eq. (24), and substituting it into Eq. (55) gives Eq. (25). The two limits follow by inspection. At 𝑑𝑑𝑥𝑥 = 0 the exponential equals one and 1 − 2𝑝𝑝 = 0, so 𝔼𝔼(|𝑇𝑇|) = 2𝜎𝜎/√𝜋𝜋 and Var�𝑋𝑋(1) � = 𝜎𝜎 2 (1 − 1/𝜋𝜋); as 𝑑𝑑𝑥𝑥 → ∞ the exponential vanishes and 1 − 2𝑝𝑝 → 1, so 𝔼𝔼(|𝑇𝑇|) → 𝑑𝑑𝑥𝑥 and Var�𝑋𝑋(1) � → 𝜎𝜎 2. Moreover Var�𝑋𝑋(1) � is nondecreasing 1 in 𝑑𝑑𝑥𝑥 , since d 𝔼𝔼(|𝑇𝑇|)/d 𝑑𝑑𝑥𝑥 = 1 − 2𝑝𝑝 and hence d Var�𝑋𝑋(1) �/d 𝑑𝑑𝑥𝑥 = [𝑑𝑑𝑥𝑥 − 𝔼𝔼(|𝑇𝑇|) (1 − 2𝑝𝑝)] =

1 2

Cov�|𝑇𝑇|, sign(𝑇𝑇)� ≥ 0, where the covariance is nonnegative because 𝔼𝔼[sign(𝑇𝑇) ∣

|𝑇𝑇| = 𝑟𝑟] = tanh(𝑑𝑑𝑥𝑥 𝑟𝑟/𝜏𝜏 2 ) is nondecreasing in 𝑟𝑟. Hence 𝜎𝜎 2 (1 − 1/𝜋𝜋) ≤ Var�𝑋𝑋(1) � ≤ 𝜎𝜎 2. █ A.4 Proof of Corollary 2

By Eqs. (14) and (15) the estimator 𝜽𝜽ˆ B minimizes the inner integral of Eq. (14) with 𝐿𝐿 = 𝐿𝐿M for every 𝑉𝑉, so it minimizes 𝑅𝑅 under 𝐿𝐿M over all functions of the frame. Both 𝜽𝜽ˆ S and 𝜽𝜽ˆ MMSE are functions of the frame. █ Funding

Acknowledgments Disclosures

The authors declare no conflicts of interest. Data Availability

Data underlying the results presented in this paper, i.e., the paper draft and all the Python scripts built on SMLM_Lib that generate the datasets, train the networks and produce the tables, are available in Ref. [39]. The SMLM_Lib library is available in Ref. [38]. References 1. 2. 3. 4.

5. 6. 7. 8. 9.

E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” Science, 313(5793), 1642–1645(2006). M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM),” Nat. Methods, 3(10), 793-796(2006). S. T. Hess, T. P. Girirajan, and M. D. Mason, “Ultra-high resolution imaging by fluorescence photoactivation localization microscopy,” Biophys. J., 91(11), 4258–4272(2006). M. Heilemann , S. v. Linde , M. Schüttpelz , R. Kasper, B. Seefeldt, A. Mukherjee, P. Tinnefeld, and M. Sauer, “Subdiffraction‐resolution fluorescence imaging with conventional fluorescent probes,” Angew. Chem. Int. Ed., 47(33), 6172–6176(2008). R. J. Ober, S. Ram and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J., 86(2), 1185-1200(2004). I. K. Mortensen, L. S. Churchman, J. A. Spudich and H. Flyvbjerg, “Optimized localization analysis for singlemolecule tracking and super-resolution microscopy,” Nat. Methods, 7(5), 377-381(2010). B. Rieger and S. Stallinga, “The lateral and axial localization uncertainty in super‐resolution light microscopy,” ChemPhysChem, 15(4), 664-670(2014). S. Stallinga and B. Rieger, “Accuracy of the Gaussian point spread function model in 2D localization microscopy,” Opt. Express, 18(24), 24461-24476(2010). B. Huang, W. Wang, M. Bates and X. Zhuang, “Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy,” Science, 319(5864), 810-813(2008).

10. S. R. P. Pavani, M. A. Thompson, J. S. Biteen, S. J. Lord, N. Liu, R. J. Twieg, R. Piestun and W. E. Moerner, “Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function,” Proc. Natl. Acad. Sci. USA, 106(9), 2995-2999(2009). 11. Y. Shechtman, S. J. Sahl and A. S. Backer, “Optimal point spread function design for 3D imaging,” Phys. Rev. Lett., 113(13), 133902(2014). 12. Y. Li, M. Mund, P. Hoess, J. Deschamps, U. Matti, B. Nijmeijer, V. J. Sabinina, J. Ellenberg, I. Schoen and J. Ries, “Real-time 3D single-molecule localization using experimental point spread functions," Nat. Methods, 15(5), 367-369(2018). 13. Y. Sun, “Information sufficient segmentation and signal-to-noise ratio in stochastic optical localization nanoscopy,” Opt. Lett., 45(21), 6102-6105(2020). 14. M. Sun and Y. Sun, “Information sufficient segmentation and signal‐to‐noise ratio for 3D astigmatism stochastic optical localization nanoscopy,” Electron. Lett., 58(2), 58-60(2022). 15. Y. Sun, “Localization precision of stochastic optical localization nanoscopy using single frames,” J. Biomed. Optics, 18(11), 111418-14(2013). 16. J. Chao, E. S. Ward and R. J. Ober, “Fisher information theory for parameter estimation in single molecule microscopy: tutorial,” JOSA A, 33(7), B36-B57(2016). 17. Y. Sun, and Y. Guan, “Effect of unknown emitter intensities on localization accuracy in stochastic optical localization nanoscopy using single frames,” JOSA A, 38(12), 1830-1840(2021). 18. Y. Sun, “Root mean square minimum distance as a quality metric for stochastic optical localization nanoscopy images,” Sci. Reports, 8(1), 17211(2018). 19. Y. Sun, “Asymptotic efficiency of the global maximum-likelihood estimator in multi-emitter localization microscopy,” arXiv: 2607.28985(2026). 20. A. P. Dempster, N. M. Laird and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. R. Stat. Soc. Ser. B, 39(1), 1-22(1977). 21. Y. Sun, Y. Zeng, W. Yen, J. M. Tarbell and B. M. Fu, “Expectation maximization can outperform Cramer-Rao lower bound in stochastic optical localization nanoscopy,” Quantitative Bioimaging Conf. (QBI2014) (2014). 22. E. L. Lehmann and G. Casella, Theory of point estimation, New York, NY: Springer (1998). 23. E. Nehme, L. E. Weiss, T. Michaeli and Y. Shechtman, “Deep-STORM: super-resolution single-molecule microscopy by deep learning,” Optica, 5(4), 458-464(2018). 24. A. Speiser, L.-R. Müller, P. Hoess, U. Matti, C. J. Obara, W. R. Legant, A. Kreshuk, J. H. Macke, J. Ries and S. C. Turaga, “Deep learning enables fast and dense single-molecule localization with high accuracy,” Nat. Methods, 18(9), 1082-1090(2021). 25. P. Zhang, S. Liu, A. Chaurasia, D. Ma, M. J. Mlodzianoski, E. Culurciello and F. Huang, “Analyzing complex single-molecule emission patterns with deep learning,” Nat. Methods, 15(11), 913-916(2018). 26. E. Nehme, D. Freedman, R. Gordon, B. Ferdman, L. E. Weiss, O. Alalouf, T. Naor, R. Orange, T. Michaeli and Y. Shechtman, “DeepSTORM3D: dense 3D localization microscopy and PSF design by deep learning,” Nat. Methods, 17(7), 734-740(2020). 27. Z. Zhou, J. Wu, Z. Wang and Z.-L. Huang, “Deep learning using a residual deconvolutional network enables real-time high-density single-molecule localization microscopy,” Biomed. Opt. Express, 14(4), 18331847(2023). 28. F. Huang, T. M. Hartwich, F. E. Rivera-Molina, Y. Lin, W. C. Duim, J. J. Long, P. D. Uchil, J. R. Myers, M. A. Baird, W. Mothes, M. W. Davidson, D. Toomre and J. Bewersdorf, “Video-rate nanoscopy using sCMOS camera–specific single-molecule localization algorithms,” Nat. Methods, 10(7), 653-658(2013). 29. H. W. Kuhn, “The Hungarian method for the assignment problem,” Nav. Res. Logist. Q., 2(1-2), 83-97(1955). 30. Y. Sun, “Partition of estimated locations: an approach to accurate quality metrics for stochastic optical localization Nanoscopy,” J. Opt. Soc. Am. A, 39(12), 2307-2315(2022). 31. F. Agostinelli, M. Hoffman, P. Sadowski and P. Baldi, “Learning activation functions to improve deep neural networks,” arXiv:1412.6830(2014). 32. S. Scardapane, S. V. Vaerenbergh, S. Totaro and A. Uncini, “Kafnets: Kernel-based non-parametric activation functions for neural networks,” Neural Networks , 110, 19-32(2019). 33. Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljačić, T. Y. Hou and M. Tegmark, “KAN: Kolmogorov-Arnold Networks,” arXiv:2404.19756(2024). 34. Y. LeCun, L. Bottou, Y. Bengio and H. Patrick, “Gradient-based learning applied to document recognition,” Proc. IEEE, 86(11), 2278-2324(1998). 35. A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit and N. Houlsby, “An image is worth 16x16 words: Transformers for image recognition at scale,” arXiv:2010.11929(2020). 36. D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv:1412.6980(2014). 37. I. Loshchilov and F. Hutter, “SGDR: Stochastic gradient descent with warm restarts,” arXiv:1608.03983(2016). 38. SunCCNY, SMLM_Lib, [Online]. Available: https://github.com/SunCCNY/Nanoscopy/tree/main/SMLM_Lib. [Accessed 5 July 2026]. 39. SunCCNY, NN-Bayes. [Online]. Available: https://github.com/SunCCNY/Nanoscopy/tree/main/2026-NNBayes. [Accessed 15 Sept. 2026]. 40. H. Ma and Y. Liu, “Super-resolution localization microscopy: Toward high throughput, high quality, and low cost,” APL photonics, 5(6), 060902(2020).

Record · ID 978441 · SHA-256 ba104150dd4ebeb4
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.