Simon Labouesse

Marc Allain

Member, IEEE Jérôme Idier

Awoke Negash

Thomas Mangeat
email: thomas.mangeat@univ-tlse3.fr

Penghuan Liu

Anne Sentenac

Sébastien Bourguignon

Joint reconstruction strategy

Keywords: Super-resolution, fluorescence microscopy, speckle imaging, near-black object model, proximal splitting

come

I. INTRODUCTION

In Structured Illumination Microscopy (SIM), the sample, characterized by its fluorescence density ρ, is illuminated successively by M distinct inhomogeneous illuminations Im. Fluorescence light emitted by the sample is collected by a microscope objective and recorded on a camera to form an image ym. In the linear regime, and with a high photon counting rate 1 , the dataset {ym} M m=1 is related to the sample ρ via [START_REF] Goodman | Introduction to Fourier Optics[END_REF]

ym = H ⊗ (ρ × Im) + εm, m = 1 • • • M, ( 1
)
where ⊗ is the convolution operator, H is the microscope point spread function (PSF) and εm is a perturbation term accounting for (electronic) noise in the detection and modeling errors. Since the spatial spectrum of the PSF [i.e., the optical transfer function (OTF)] is strictly bounded by its cut-off frequency, say, νpsf, if the illumination pattern Im is homogeneous, then the spatial spectrum of ρ that can be retrieved from the image ym is restricted to frequencies below νpsf. When the illuminations are inhomogeneous, frequencies beyond νpsf can be recovered from the low resolution images because the illuminations, acting as carrier waves, downshift part of the spectrum inside the OTF support [START_REF] Heintzmann | Laterally modulated excitation microscopy: improvement of resolution by using a diffraction grating[END_REF], [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF]. Standard SIM resorts to harmonic illumination patterns for which the reconstruction of the super-resolved image can be easily done by solving a linear system in the Fourier domain. In this case, the gain in resolution depends on the OTF support, the illumination cut-off frequency and the available signal-to-noise ratio (SNR). The main drawback of SIM is that it requires the knowledge of the illumination patterns and thus a stringent control of the experimental setup. If these patterns are not known with sufficient accuracy [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF], [START_REF] Ayuk | Structured illumination fluorescence microscopy with distorted excitations using a filtered blind-SIM algorithm[END_REF], severe artifacts appear in the reconstruction. Specific estimation techniques have been developed for retrieving the parameters of the periodic patterns from the images [START_REF] Orieux | Bayesian estimation for optimized structured illumination microscopy[END_REF]- [START_REF] Wicker | Non-iterative determination of pattern phase in structured illumination microscopy using auto-correlations in Fourier space[END_REF], but they can fail if the SNR is too low or if the excitation patterns are distorted, e.g., by inhomogeneities in the sample refraction index. The Blind-SIM strategy [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF], [START_REF] Ayuk | Structured illumination fluorescence microscopy with distorted excitations using a filtered blind-SIM algorithm[END_REF], [START_REF] Negash | Improving the axial and lateral resolution of three-dimensional fluorescence microscopy using random speckle illuminations[END_REF] has been proposed to tackle this key issue, the principle being to retrieve the sample fluorescence density without the knowledge of the illumination patterns. In addition, speckle illumination patterns are promoted instead of harmonic ones, the latter being much more difficult to generate and control. From the methodological viewpoint, this strategy relies on the simultaneous (joint) reconstruction of the fluorescence density and of the illumination patterns. More precisely, joint reconstruction is achieved through the iterative resolution of a constrained least-squares problem. However, the computational time of such a scheme clearly restricts the applicability of the method.

This paper provides a global re-foundation of the joint Blind-SIM strategy. More specifically, our work develops two specific, yet complementary, contributions:

• The joint Blind-SIM reconstruction problem is first revisited, resulting in an improved numerical implementation with execution times decreased by several orders of magnitude. Such an acceleration relies on two technical contributions. Firstly, we show that the problem proposed in [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF] is equivalent to a fully separable constrained minimization problem, hence bringing the original (large-scale) problem to M sub-problems with smaller scales. Then, we introduce a new preconditioned proximal iteration (denoted PPDS) to efficiently solve each sub-problem. The PPDS strategy is an important contribution of this article: it is provably convergent [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF], easy to implement and, for our specific problem, we empirically observe a superlinear asymptotic convergence rate. With these elements, the joint Blind-SIM reconstruction proposed in this paper is fast and can be highly parallelized, opening the way to real-time reconstructions.

• Beside these algorithmic issues, the mechanism driving superresolution (SR) in this blind context is investigated, and a connection is established with the well-known "Near-black object" effect introduced in Donoho's seminal contribution [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF]. We show that the SR relies on sparsity and positivity constraints enforced by the unknown illumination patterns. This finding helps to understand in which situation super-resolved reconstructions may be provided or not. A significant part of this work is then dedicated to numerical simulations aiming at illustrating how the SR effect can be enhanced. In this perspective, our simulations show that two-photon speckle illuminations potentially increase the SR power of the proposed method.

The pivotal role played by sparse illuminations in this SR mechanism also draws a connexion between joint Blind-SIM and other random activation strategies like PALM [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF] or STORM [START_REF] Rust | Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM)[END_REF]; see also [START_REF] Mukamel | Statistical deconvolution for superresolution fluorescence microscopy[END_REF], [START_REF] Min | FALCON: fast and unbiased reconstruction of high-density super-resolution microscopy data[END_REF] for explicit sparsity methods applied to STORM. With PALM/STORM, unparalleled resolutions result from an activation process that is massively sparse and mostly localized on the marked structures. With the joint Blind-SIM strategy, the illumination pattern playing the role of the activation process is not that "efficient" and lower resolutions are obviously expected. Joint Blind-SIM however provides SR as long as the illumination patterns enforce many zero (or almost zero) values in the product ρ × Im: the sparser the illuminations, the higher the expected resolution gain with joint Blind-SIM. Such super resolution can be induced by either deterministic or random patterns. Let us mention that random illuminations are easy and cheap to generate, and that a few recent contributions advocate the use of speckle illuminations for super-resolved imaging, either in fluorescence [START_REF] Min | Fluorescent microscopy beyond diffraction limits using speckle illumination and joint support recovery[END_REF], [START_REF] Oh | Sub-Rayleigh imaging via speckle illumination[END_REF] or in photo-acoustic [START_REF] Chaigne | Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination[END_REF] microscopy. In these contributions, however, the reconstruction strategies are derived from the statistical modeling of the speckle, hence, relying on the random character of the illumination patterns. In comparison, our approach only requires that the illuminations cancel-out the fluorescent object and that their sum is known with sufficient accuracy. Finally, we also note that [START_REF] Labouesse | Fluorescence blind structured illumination microscopy: a new reconstruction strategy[END_REF] corresponds to an early version of this work. Compared to [START_REF] Labouesse | Fluorescence blind structured illumination microscopy: a new reconstruction strategy[END_REF], several important contributions are presented here, mainly: the super-resolving power of Blind-SIM is now studied in details, and a comprehensive presentation of the proposed PPDS algorithm includes a tuning strategy for the algorithm parameter that allows a substantial reduction of the computation time.

The remainder of the paper is organized as follows. In Section II, the original Blind-SIM formulation is introduced and further simplified; this reformulation is then used to get some insight on the mechanism that drives the SR in the method. Taking advantage of this analysis, a penalized Blind-SIM strategy is proposed and characterized with synthetic data in Section III. Finally, the PPDS algorithm developed to cope with the minimization problem is presented and tested in Section IV, and conclusions are drawn in Section V.

II. SUPER-RESOLUTION WITH JOINT BLIND-SIM ESTIMATION

In the sequel, we focus on a discretized formulation of the observation model [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF]. Solving the two-dimensional (2D) Blind-SIM reconstruction problem is equivalent to finding a joint solution ( ρ, { Im} M m=1 ) to the following constrained minimization problem [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF]:

min ρ,{Im} M m=1 ym -Hdiag(ρ) Im 2 (2a) subject to m Im = M × I0 (2b) and ρn ≥ 0, Im;n ≥ 0, ∀m, n (2c)
with H ∈ R P ×N the 2D convolution matrix built from the discretized PSF. We also denote ρ = vect(ρn) ∈ R N the discretized fluorescence density, ym = vect(ym;n) ∈ R P the m-th recorded image, and Im = vect(Im;n) ∈ R N the m-th illumination with expected spatial intensity I0 = vect(I0;n) ∈ R N + (this latter quantity may be spatially inhomogeneous but it is supposed to be known). Let us remark that (2) is a biquadratic problem. Block coordinate descent alternating between the object and the illuminations could be a possible minimization strategy, relying on cyclically solving M + 1 quadratic programming problems [START_REF] Jost | Optical sectioning and high resolution in single-slice structured illumination microscopy by thick slice blind-SIM reconstruction[END_REF]. In [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF], a more efficient but more complex scheme is proposed. However, the minimization problem (2) has a very specific structure, yielding a fast and simple strategy, as shown below.

A. Reformulation of the optimization problem

According to [START_REF] Labouesse | Fluorescence blind structured illumination microscopy: a new reconstruction strategy[END_REF], let us first consider problem (2) without the equality constraint (2b). It is equivalent to M independent quadratic minimization problems:

minq m ym -Hqm 2 (3a) subject to qm ≥ 0, (3b)
where we set qm := vect(ρn × Im;n). Each minimization problem (3) can be solved in a simple and efficient way (see Sec. IV), hence providing a set of global minimizers { qm} M m=1 . Although the latter set corresponds to an infinite number of solutions ( ρ, { Im} M m=1 ), the equality constraint (2b) defines a unique solution such that qm = vect( ρn × Im;n) for all m:

ρ = Diag(I0) -1 q (4a) ∀m Im = Diag( ρ) -1 qm (4b)
with q := 1 M m qm. The solution (4) exists as long as I0;n = 0 and ρn = 0, ∀n. The first condition is met if the sample is illuminated everywhere (in average), which is an obvious minimal requirement. For any pixel sample such that ρn = 0, the corresponding illumination Im;n is not defined; this is not a problem as long as the fluorescence density ρ is the only quantity of interest. Let us also note that the following implication holds:

I0;n ≥ 0, qm;n ≥ 0 =⇒ Im,n ≥ 0 and ρn ≥ 0.
Because we are dealing with intensity patterns, the condition I0;n ≥ 0 is always met, hence the positivity granted for both the density and the illumination estimates, i.e., the positivity constraint (2c), is granted by [START_REF] Heintzmann | Laterally modulated excitation microscopy: improvement of resolution by using a diffraction grating[END_REF]. Indeed, it should be clear that combining [START_REF] Goodman | Introduction to Fourier Optics[END_REF] and (4) solves the original minimization problem (2): on the one hand, the equality constraint (2b) is met since2

m Im = Diag( ρ) -1 q = M I0 (5)
and on the other hand, the solution (4) minimizes the criterion given in (2a) since it is built from { qm} M m=1 , which minimizes (3a). Finally, it is worth noting that the constrained minimization problem (2) may have multiple solutions. In our reformulation, this ambiguity issue arises in the "minimization step" (3): while each problem (3) is convex quadratic, and thus admits only global solutions (which in turn provide a global solution to problem (2) when recombined according to (4a)-(4b)), it may not admit unique solutions since each criterion (3a) is not strictly convex3 in qm. Furthermore, the positivity constraint (3b) prevents any direct analysis of these ambiguities. The next subsection underlines however the central role of this constraint in the joint Blind-SIM strategy originally proposed in [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF].

B. Super-resolution unveiled

Whereas the mechanism that conveys SR with known structured illuminations is well understood (see [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF] for instance), the SR capacity of joint blind-SIM has not been characterized yet. It can be made clear, however, that the positivity constraint (2c) plays a central role in this regard. Let H + be the pseudo-inverse of H [START_REF] Golub | Matrix computation[END_REF]Sec. 5.5.4]. Then, any solution to the problem (2a)-(2b), i.e, without positivity constraints, reads

ρ = Diag(I0) -1 (H + y + q ⊥ ) (6a) Im = Diag( ρ) -1 (H + ym + q ⊥ m ), (6b)
with y = 1 M m ym, and q ⊥ = 1 M m q ⊥ m where q ⊥ m is an arbitrary element of the kernel of H, i.e. with arbitrary frequency components above the OTF cutoff frequency. Hence, the formulation (2a)-(2b) has no capacity to discriminate the correct high frequency components, which means that it has no SR capacity. Under the positivity constraint (2c), we thus expect that the SR mechanism rests on the fact that each illumination pattern Im activates the positivity constraint on qm in a frequent manner.

A numerical experiment is now considered to support this assertion. A set of M collected images are simulated following [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF] with the PSF H given by the usual Airy pattern that reads in polar coordinates

H(r, θ) = k 2 0 π J1(r k0 NA) k0 r 2 , r ≥ 0, θ ∈ R, ( 7
)
where J1 is the first order Bessel function of the first kind, NA is the objective numerical aperture set to 1.49, and k0 = 2π/λ is the free-space wavenumber with λ the emission/excitation wavelength. The ground truth is the 2D 'star-like' fluorescence pattern depicted in Fig. 1(left). The image sampling step for all the simulations involving the star pattern is set 4 to λ/20. For this numerical simulation, the illumination set {Im} M m=1 consists in M = 200 modified speckle patterns, see Fig. 2(A). More precisely, a first set of illuminations is 4 For an optical system modeled by [START_REF] Orieux | Bayesian estimation for optimized structured illumination microscopy[END_REF], the sampling rate of the (diffraction-limited) acquisition is usually the Nyquist rate driven by the OTF cutoff frequency ν psf = 2k 0 NA. A higher sampling rate is obviously needed for the super-resolved reconstruction, the up-sampling factor between the "acquisition" and the "processing" rates being at least equal to the expected SR factor. Here, we adopt a common sampling rate for any simulation involving the star-like pattern (even with diffraction-limited images), as it allows a direct comparison of the reconstruction results. obtained by adding a positive constant (equal to 3) to each speckle pattern, resulting in illuminations that never activate the positivity constraint in (3). On the contrary, the second set of illuminations is built by subtracting a small positive constant (equal to 0.2) to each speckle pattern, the negative values being set to zero. The resulting illuminations are thus expected to activate the positivity constraint in [START_REF] Goodman | Introduction to Fourier Optics[END_REF]. For both illumination sets, low-resolution microscope images are simulated and corrupted with Gaussian noise; in this case, the standard deviation was chosen so that the SNR of the total dataset is 40 dB. Corresponding reconstructions of the first product image q1 obtained via the resolution of ( 3) is shown in Fig. 2(B), while the retrieved sample (4a) is shown in Fig. 2(C); for each reconstruction, the spatial mean I0 in (4a) is set to the statistical expectation of the corresponding illumination set. As expected, the reconstruction with the first illumination set is almost identical to the deconvolution of the wide-field image shown in Fig. 1(upper-right), i.e., there is no SR in this case. On the contrary, the second set of illuminations produces a super-resolved reconstruction, hence establishing the central role of the positivity constraint in the original joint reconstruction problem (2).

III. A PENALIZED APPROACH FOR JOINT BLIND-SIM

As underlined in the beginning of Subsection II-B, there is an ambiguity issue concerning the original joint Blind-SIM reconstruction problem. A simple way to enforce unicity is to slightly modify (3) by adding a strictly convex penalization term. We are thus led to solving

min qm≥0 ym -Hqm 2 + ϕ(qm). ( 8
)
Another advantage of such an approach is that ϕ can be chosen so that robustness to the noise is granted and/or some expected features in the solution are enforced. In particular, the analysis conveyed above suggests that favoring sparsity in each qm is suited since speckle or periodic illumination patterns tend to frequently cancel or nearly cancel the product images qm. For such illuminations, the Near-Black Object introduced in Donoho's seminal paper [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF] is an appropriate modeling and, following this line, we found that the separable " 1 +

2" penalty 5 provides super-resolved reconstructions:

ϕ(qm) := α n |qm;n| + β||qm|| 2 , α ≥ 0, β > 0. ( 9
)
With properly tuned (α, β), our joint Blind-SIM strategy is expected to bring SR if "sparse" illumination patterns Im are used, i.e., if they enforce qm;n = 0 for most (or at least many) n. More specifically, it is shown in [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF]Sec. 4] that SR occurs if the number of non-zero Im;n (i.e., the number of non-zero components to retrieve in qm) divided by N is lower than 1 2 R/N , with R/N the incompleteness ratio and R the rank of H. In addition, the resolving power is driven by the spacing between the components to retrieve that, ideally, should be greater than the Rayleigh distance λ 2 NA , see [12, pp. 56-57]. These conditions are rather stringent and hardly met by illumination patterns that can be reasonably considered in practice. These illumination patterns are usually either deterministic harmonic or quasi-harmonic 6 patterns, or random speckle patterns, these latter illuminations being much easier to generate [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF]. Nevertheless, in both cases, a SR effect is observed in joint Blind-SIM. Moreover, one can try to maximize this effect via the tuning of some experimental parameters that are left to the designer of the setup. Such parameters are mainly: the period of the light grid and the number of grid shifts for harmonic patterns, the spatial correlation length and the point-wise statistics of the speckle patterns. Investigating the SR properties with respect to these parameters on a theoretical ground seems out of reach. However, a numerical analysis is possible and some illustrative results are now provided that address this question. Reconstructions shown in the sequel are built from (4a) via the numerical resolution of ( 8)- [START_REF] Wicker | Non-iterative determination of pattern phase in structured illumination microscopy using auto-correlations in Fourier space[END_REF]. For sake of clarity, all the algorithmic details concerning this minimization problem are reported in Sec. IV. These simulations were performed with low-resolution microscope images corrupted by additive Gaussian noise such that the signal-to-noise ratio (SNR) of the dataset {ym} M m=1 is 40 dB. In addition, we note that this penalized joint Blind-SIM strategy requires an explicit tuning of some hyper-parameters, namely α and β in the regularization function [START_REF] Wicker | Non-iterative determination of pattern phase in structured illumination microscopy using auto-correlations in Fourier space[END_REF]. Further details concerning these parameters are reported in Sec. III-D.

A. Regular and distorted harmonic patterns

We first consider unknown harmonic patterns defining a "standard" SIM experiment with M = 18 patterns. More precisely, the illuminations are harmonic patterns of the form I(r) = 1 + cos(2πν t r + φ) where φ is the phase shift, and with r = (x, y) t and ν = (νx, νy) t the spatial coordinates and the spatial frequencies of the harmonic 5 The super-resolved solution in [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF] is obtained with a positivity constraint and a 1 separable penalty. However, ambiguous solutions may exist in this case since the criterion to minimize is not strictly convex. The 2 penalty in ( 9) is then mostly introduced for the technical reason that a unique solution exists for problem [START_REF] Wicker | Phase optimisation for structured illumination microscopy[END_REF]. 6 Dealing with distorted patterns is of particular practical importance since it allows to cope with the distortions and misalignments induced by the instrumental uncertainties or even by the sample itself [START_REF] Ayuk | Structured illumination fluorescence microscopy with distorted excitations using a filtered blind-SIM algorithm[END_REF], [START_REF] Jost | Optical sectioning and high resolution in single-slice structured illumination microscopy by thick slice blind-SIM reconstruction[END_REF]. function, respectively. Distorted versions of these patterns (deformed by optical aberrations such as astigmatism and coma) were also considered. Three distinct orientations θ := tan -1 (νy/νx) ∈ {0, 2π/3, 4π/3}, for each of which six phase shifts of one sixth of the period, were considered. The frequency of the harmonic patterns ||ν|| := (ν 2 x + ν 2 y ) 1/2 is set to 80% of the OTF cutoff frequency, i.e., it lies inside the OTF support. One regular and one distorted pattern are depicted in Fig. 3(A) and the penalized joint Blind-SIM reconstructions are shown in Fig. 3(B). For both illumination sets, a clear SR effect occurs, which is similar to the one obtained with the original approach presented in [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF]. As expected, however, the reconstruction quality achieved in this blind context is lower than what can be obtained with standard harmonic SIM -for the sake of comparison, see Fig. 1(B). In addition, we note that some artifacts may appear if the number of phase shifts for each orientation is decreased, see Fig. 3(C-left). If we keep in mind that the retrieved

B. Speckle illumination patterns

We now consider second-order stationary speckle illuminations Im with known first order statistics I0;n = I0, ∀n. Each one of these patterns is a fully-developed speckle drawn from the pointwise intensity of a correlated circular Gaussian random field. The correlation is adjusted so that the pattern Im exhibits a spatial correlation of the form (7) but with "numerical aperture" parameter NAill that sets the correlation length to λ 2 NA ill within the random field. As an illustration, the speckle pattern shown in Fig. 4(A-left) was generated in the standard case 7 NAill = NA. From this set of regular (fully-developed) speckle patterns, we also consider another set of random illumination patterns built by squaring each speckle pattern, see Fig. 4(A). These "squared" patterns are considered hereafter because they give a deeper insight about the SR mechanism at work in joint Blind-SIM. Moreover, we discuss later in this subsection that 7 It is usually considered that NA ill = NA if the illumination and the collection of the fluorescent light are performed through the same optical device. these patterns can be generated with other microscopy techniques, hence extending the concept of random illumination microscope to other optical configurations. From a statistical viewpoint, the probability distribution function (pdf) of "standard" and "squared" speckle patterns differ. For instance, the pdf of the squared speckle intensity is more concentrated around zero8 than the exponential pdf of the standard speckle intensity. In addition, the spatial correlation is also changed since the power spectral density of the "squared" random field spans twice the initial support of its speckle counterpart [START_REF] Denk | Two-photon laser scanning fluorescence microscopy[END_REF]. As a result, the "squared" speckle grains are sharper, and they enjoy larger spatial separation. According to previous SR theoretical results [12, p. 57] (see also the beginning of Sec. III), these features may bring more SR in joint Blind-SIM than standard speckle patterns. This assumption was indeed corroborated by our simulations. For instance, the reconstructions in Fig. 4(B) were obtained from a single set of M = 1000 speckle patterns such that NAill = NA: in this case, the "squared" illuminations (obtained by squaring the speckle patterns) provide a higher level of SR than the standard speckle illuminations.

Figure 5 shows how the reconstruction quality varies with the number of illumination patterns. With very few illuminations, the sample is retrieved in the few places that are activated by the "hot spots" of the speckle patterns. This actually illustrates that the joint Blind-SIM approach is also an "activation" strategy in the spirit of PALM [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF] and STORM [START_REF] Rust | Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM)[END_REF]. With our strategy, the activation process is nevertheless enforced by the structured illumination patterns and (A) not by the fluorescent markers staining the sample. This effect is more visible with the squared illumination patterns and, with these somehow sparser illuminations, the number of patterns needs to be increased so that the fluctuations in m Im is moderate, hence making the equality (2b) a legitimate constraint. We also stress that these simulations corroborate the empirical statement that M ≈ 9 harmonic illuminations and M ≈ 200 speckle illuminations produce comparable super-resolved reconstructions, see Fig. 3(Cleft) and Fig. 5(B-left). Obviously, imaging with random speckle patterns remains an attractive strategy since it is achieved with a very simple experimental setup, see [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF] for details. For both random patterns, we also note that increasing the correlation length above the Rayleigh distance λ 2 NA (i.e., setting NAill < NA) deteriorates the SR whereas, conversely, taking NAill = 2NA enhances it, see Fig. 6-(A,B). However, the resolving power of the joint Blind-SIM estimate deteriorates if the correlation length is further decreased; for instance, uncorrelated speckle patterns are finally found to hardly produce any SR, see Fig. 6-(C). Indeed, with arbitrary small correlation lengths, many "hot spots" tend to be generated within a single Rayleigh distance, leading to this loss in the resolving power. Obviously, the "squared" speckle patterns are less sensitive to this problem because they are inherently sparser.

Finally, the experimental relevance of the simulations involving "squared" speckle illuminations needs to be addressed. Since a twophoton (2P) fluorescence interaction is sensitive to the square of the intensity [START_REF] Gu | Comparison of three-dimensional imaging properties between two-photon and single-photon fluorescence microscopy[END_REF], most of these simulations can actually be considered as wide-field 2P structured illumination experiments. Unlike one-photon (i.e., fully-developed) speckle illuminations 9 , though, a 2P interaction requires an excitation wavelength λill ∼ 1000 nm that is roughly twice the one of the collected fluorescence λdet ∼ 500 nm. The lateral 2P correlation length being λ ill 4NA ill , epi-illumination setups with onephoton (1P) and 2P illuminations provide similar lateral correlation lengths. This 2P instrumental configuration is simulated in Fig. 6(Aright), which does not show any significant SR improvement with respect to 1P epi-illumination interaction shown in Fig. 5(C-left). The increased SR effect driven by "squared" illumination patterns can nevertheless be obtained with 2P interactions if the excitation and the collection are performed with separate objectives. For instance, the behaviors shown in Fig. 5(C-right) and in Fig. 6(B-right) can be obtained if the excitation NA is, respectively, twice and four times the collection NA. With these configurations, the 2P excitation exhibits a correlation length which is significantly smaller than the one driven by the objective PSF, and a strong SR improvement is observed in simulation by joint Blind-SIM. The less spectacular simulation shown in Fig. 6(C-right) can also be considered as a 2P excitation, in the "limit" case of a very low collection NA. The 1P simulation shown in Fig. 6(C-left) rather mock a photo-acoustic imaging experiment [START_REF] Chaigne | Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination[END_REF], an imaging technique for which the illumination lateral correlation length is negligible with respect to the PSF width.

As a final remark, we stress that 2P interactions are not the only way to generate sparse illumination patterns for the joint Blind-SIM. In particular, super-Rayleigh speckle patterns [START_REF] Bromberg | Generating non-Rayleigh speckles with tailored intensity statistics[END_REF] are promising candidates for that purpose.

C. Some reconstructions from real and mock data

The star test-pattern used so far is a simple and legitimate mean to evaluate the resolving power of our strategy [START_REF] Horstmeyer | Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM)[END_REF], but it hardly provides a convincing illustration of what can be expected with real data. Therefore, we now consider the processing of more realistic datasets with joint Blind-SIM. In this section, the microscope acquisitions are all designed so that the spatial sampling rate is equal or slightly above the Nyquist rate λ 4NA . As a consequence, a preliminary up-sampling step of the camera acquisitions is performed so that their sampling rate reaches that of the super-resolved reconstruction.

As a first illustration, we consider a real dataset resulting from a test sample composed of fluorescent beads with diameters of 100 nm. A set of 100 one-photon speckle patterns is generated by a spatial light modulator and a laser source operating at λill = 488 nm. The fluorescent light at λcoll = 520 nm is collected through an objective with NA = 1.49 and recorded by a camera. The excitation and the collection of the emitted light are performed through the same objective, i.e., the setup is in epi-illumination mode. The total number of photons per camera pixels is about 65 000. In the perspective of further processing, this set of camera acquisitions is first up-sampled with a factor of two. Figure 7(A-left) shows the sum of these (upsampled) acquisitions, which is similar to a wide-field image. Wiener deconvolution of this image can be performed so that all spatial frequencies transmitted by the OTF are equivalently contributing in a diffraction-limited image of the beads, see Figure 7(A-middle). The processing of the dataset by the joint Blind-SIM strategy shown in Figure 7(A-right) reveals several beads that are otherwise unresolved on the diffraction-limited images, hence demonstrating a clear SR effect. In this case, the distance between the closest pair of resolved beads provides an upper bound for the final resolution, that is λcoll/5. 9 With one-photon interactions, the Stokes shift [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF] implies that the excitation and the fluorescence wavelengths are not strictly equivalent. The difference is however negligible in practice (about 10%), hence our assumption that one-photon interactions occur with identical wavelengths for both the excitation and the collection. The experimental demonstration above does not involve any biological sample, and we now consider a simulation designed to be close to a real-world biological experiment. More specifically, the STORM reconstruction of a marked neuron 10 is used as a ground truth to simulate a series of microscope acquisitions generated from onephoton speckle illuminations. Our simulation considers 300 illuminations and acquisitions, both performed through the same objective, at λ = 488 nm and with NA = 1. Each low-resolution acquisition is finally plagued with Poisson noise, the total photon budget being equal to 50 000 so that it fits to the one of a standard fluorescence wide-field image. The sample (ground truth) shown in Figure 7(Bleft) interestingly exhibits a lattice structure with a 190 nm periodicity (in average) that is not resolved by the diffracted-limited image shown in Figure 7(B-middle). The joint Blind-SIM reconstruction in Figure 7(A-right) shows a significant improvement of the resolution, which reveals some parts of the underlying structure.

D. Tuning the regularization parameters

The tuning of parameters α and β in ( 9) is a pivotal issue since inappropriate values result in deteriorated reconstructions. On the one hand, the quadratic penalty in [START_REF] Wicker | Non-iterative determination of pattern phase in structured illumination microscopy using auto-correlations in Fourier space[END_REF] was mostly introduced to ensure that the minimizer defined by ( 8) is unique (via strict convexity of the criterion). However, because high-frequency components in qm are progressively damped as β increases, the latter parameter can 10 A rat hippocampal neuron in culture labelled with an anti-βIV-spectrin primary and a donkey anti-rabbit Alexa Fluor 647 secondary antibodies, imaged by STORM and processed similarly to [START_REF] Leterrier | Nanoscale architecture of the axon initial segment reveals an organized and robust scaffold[END_REF]. also be adjusted in order to prevent an over-amplification of the instrumental noise. A trade-off should nevertheless be sought since large values of β prevent super-resolution to occur. For a given SNR, β is then maintained to a fixed (usually small) value. For instance, we chose β = 10 -6 for all the simulations involving the star pattern in this paper since they were performed with a rather favorable SNR. On the other hand, the quality of reconstruction crucially depends on parameter α. More precisely, larger values of α will provide sparser solutions qm, and thus a sparser reconstructed object ρ. Fig. 8 shows an example of under-regularized and overregularized solutions, respectively corresponding to a too small and a too large value of α. The prediction of the appropriate level of sparsity to seek for each qm, or equivalently the tuning of the regularization parameter α, is not an easy task. Two main approaches can be considered. One relies on automatic tuning. For instance, a simple method called Morozov's discrepancy principle considers that the least-squares terms ym -H qm 2 should be chosen in proportion with the variance of the additive noise, the latter being assumed known [START_REF] Morozov | Methods for Solving Incorrectly Posed Problems[END_REF]. Other possibilities seek a trade-off between ym -H qm 2 and ϕ( qm). This is the case with the L-curve [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the L-curve[END_REF], but also with the recent contribution [START_REF] Song | Regularization parameter estimation for non-negative hyperspectral image deconvolution[END_REF], which deals with a situation comparable to ours. Another option relies on a Bayesian interpretation of qm as a maximum a posteriori solution, which opens the way to the estimation of α marginally of qm. In this setting, Markov Chain Monte Carlo sampling [START_REF] Lucka | Fast Markov chain Monte Carlo sampling for sparse Bayesian inference in high-dimensional inverse problems using L1-type priors[END_REF] or variational Bayes methods [START_REF] Babacan | Variational Bayesian blind deconvolution using a total variation prior[END_REF] could be employed. An alternate approach to automatic tuning consists in relying on a calibration step. It amounts to consider that similar acquisition conditions, applied to a given type of biological samples, lead to similar ranges of admissible values for the tuning of α. The validation of such a principle is however outside the scope of this article as it requires various experimental acquisitions from biological samples with known structures (or, at least, with some calibrated test patterns). Concerning the examples proposed in the present section, the much simpler strategy consisted in selecting the reconstruction which is visually the "best" among the reconstructed images with varying α.

IV. A NEW PRECONDITIONED PROXIMAL ITERATION

We now consider the algorithmic issues involved in the constrained optimization problem ( 8)-( 9). For sake of simplicity, the subscript m in ym and qm will be dropped. The reader should however keep in mind that the algorithms presented below only aim at solving one of the M sub-problems involved in the final joint Blind-SIM reconstruction. Moreover, we stress that all simulations presented in this article are performed with a convolution matrix H with a blockcirculant with circulant-block (BCCB) structure. The more general case of block-Toeplitz with Toeplitz-block (BTTB) structure is shortly addressed at the end of Subsection IV-C.

At first, let us note that ( 8)-( 9) is an instance of the more general problem min

q∈R N [f (q) := g(q) + h(q)] (10)
where g and h are closed-convex functions that may not share the same regularity assumptions: g is supposed to be a smooth function with a L-Lipschitz continuous gradient ∇g, but h does not need to be smooth. Such a splitting aims at solving constrained non-smooth optimization problems by proximal (or forward-backward) iterations.

The next subsection presents the basic proximal algorithm and the well-known FISTA that usually improves the convergence speed.

A. Basic proximal and FISTA iterations

We first present the resolution of (10) in a general setting, then the penalized joint Blind-SIM problem (8) is addressed as our particular case of interest.

1) General setting: Let q (0) be an arbitrary initial guess, the basic proximal update k → k + 1 for minimizing the convex criterion f is [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF]- [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] q (k+1) ←-P γh q (k) -γ∇g(q (k) )

where P γh is the proximity operator (or Moreau envelope) of the function γh [37, p.339]

P γh (q) := arg min x∈R N h(x) + 1 2γ ||x -q|| 2 . ( 12
)
(a) FISTA 10 it. For all these simulations, the initial-guess is q (0) = 0 and the regularization parameters is set to (α = 0.3, β = 10 -6 ).

The PPDS iteration implements the preconditioner given in [START_REF] Leterrier | Nanoscale architecture of the axon initial segment reveals an organized and robust scaffold[END_REF] with C = H t H and a = 1, see Sec IV-C for details.

Although this operator defines the update implicitly, an explicit form is actually available for many of the functions met in signal and image processing applications, see for instance [36, Table 10.2].

The Lipschitz constant L granted to ∇g plays an important role in the convergence of iterations [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]. In particular, global convergence toward a solution of (10) occurs as long as the step size γ is chosen such that 0 < γ < 2/L. However, the convergence speed is usually very low and the following accelerated version named FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] is usually preferred

q (k+1) ←-P γh ω (k) -γ∇g(ω (k) ) (13a)
ω (k+1) ←-q (k+1) + k-1 k+2 q (k+1) -q (k) . (13b)
The convergence speed toward minq f (q) achieved by ( 13) is O(1/k 2 ), which is often considered as a substantial gain compared to the O(1/k) rate of the basic proximal iteration. It should be noted however that this "accelerated" form may not always provide a faster convergence speed with respect to its standard counterpart, see for instance [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]Fig. 10.2]. FISTA was nevertheless found to be faster for solving the constrained minimization problem involved in joint Blind-SIM, see Fig. 11. We finally stress that convergence of ( 13) is granted for 0 < γ < 1/L [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

2) Solution of the m-th joint Blind-SIM sub-problem: For the penalized joint Blind-SIM problem considered in this paper, the minimization problem [START_REF] Wicker | Phase optimisation for structured illumination microscopy[END_REF] [equipped with the penalty ( 9)] takes the form [START_REF] Negash | Improving the axial and lateral resolution of three-dimensional fluorescence microscopy using random speckle illuminations[END_REF] with

g(q) = ||y -Hq|| 2 + β||q|| 2 (14a) h(q) = α n φ(qn) (14b)
where φ : R → R ∪ {+∞} is such that

φ(u) := u if u ≥ 0. +∞ otherwise. ( 14c
)
The gradient of the regular part in the splitting

∇g(q) = 2 H t (Hq -y) + βq (15)
is L-Lipschitz-continuous with L = 2 λmax(H t H) + β where λmax(A) denotes the highest eigenvalue of the matrix A. Furthermore, the proximity operator [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF] with h defined by (14b) leads to the well-known soft-thresholding rule [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized Gaussian and complexity priors[END_REF], [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF] P γh (q) = vect (max{qn -γα, 0}) .

From a practical perspective, both the basic iteration [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF] and its accelerated counterpart [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF] are easily implemented at a very low computational cost 11 from equations ( 15) and ( 16). For our penalized joint Blind-SIM approach, however, we observed that both algorithms exhibit similar convergence behavior in terms of visual aspect of the current estimate. The convergence speed is also significantly slow: several hundreds of iterations are usually required for solving the M = 200 sub-problems involved in the joint Blind-SIM reconstruction shown in Fig. 5(B). In addition, Fig. 9(ac) shows the reconstruction built with ten, fifty and one thousand FISTA iterations. Clearly, we would like that this latter quality of reconstruction is reached in a reasonable amount of time. The next subsection introduces a preconditioned primal-dual splitting strategy that achieves a much higher convergence speed, as illustrated by Fig. 9(right).

B. Preconditioned primal-dual splitting

The preconditioning technique [42, p. 69] is formally equivalent to addressing the initial minimization problem (10) via a linear transformation q := P v, where P ∈ R N ×N is a symmetric positive-definite matrix. There is no formal difficulty in defining a preconditioned version of the proximal iteration [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]. However, if one excepts the special case of diagonal matrices P [43]- [START_REF] Raguet | Preconditioning of a generalized forwardbackward splitting and application to optimization on graphs[END_REF], the proximity operator of H(v) := h(P v) cannot be obtained explicitly and needs to be computed approximately. As a result, solving a nested optimization problem is required at each iteration, hence increasing the overall computational cost of the algorithm and raising a convergence issue since the sub-iterations must be truncated in practice [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF], [START_REF] Chouzenoux | Variable metric forwardbackward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF]. Despite this difficulty, the preconditioning is widely accepted as a very effective way for accelerating proximal iterations. In the sequel, the versatile primal-dual splitting technique introduced in [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF], [START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF], [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators[END_REF] is used to propose a new preconditioned proximal iteration, without any nested optimization problem.

This new preconditioning technique is now presented for the generic problem [START_REF] Negash | Improving the axial and lateral resolution of three-dimensional fluorescence microscopy using random speckle illuminations[END_REF]. At first, we express the criterion f with respect to the transformed variables

f (P v) = G(v) + h(P v) ( 17
)
11 Since H is a convolution matrix, the computation of the gradient (15) can be performed by fast Fourier transform and vector dot-products, see for instance [START_REF] Vogel | Computational Methods for Inverse Problems[END_REF]Sec. 5.2.3].

with G(v) := g(P v). Since the criterion above is a particular case of the form considered in [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]Eq. (45)], it can be optimized by a primal-dual iteration [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]Eq. (55)] that reads

v (k+1) ←-v (k) -θτ d (k) (18a) ω (k+1) ←-ω (k) + θ∆ (k) (18b) with d (k) := ∇G(v (k) ) + P ω (k) (19a) ∆ (k) := P σh ω (k) + σP (v (k) -2τ d (k) ) -ω (k) (19b)
where the proximal mapping applied to h , the Fenchel conjugate function for h, is easily obtained from

P σh (ω) = ω -σP h/σ (ω/σ). (20)
The primal update (18a) can also be expressed with respect to the untransformed variables q:

q (k+1) ←-q (k) -θτ Bζ (k) (21)
with ζ (k) := ∇g(q (k) )+ω (k) and B := P P . Since the update ( 21) is a preconditioned primal step, we expect that a clever choice of the preconditioning matrix B will provide a significant acceleration of the primal-dual procedure. In addition, we note that the quantity

a (k) := ω (k) + σP (v (k) -2τ d (k)
) involved in the dual step via (19b) also reads

a (k) := ω (k) + σ(q (k) -2τ Bζ (k) ). (22)
Hereafter, the primal-dual updating pair (18b) and ( 21) is called a preconditioned primal-dual splitting (PPDS) iteration. Following [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]Theorem 5.1], the convergence of these PPDS iterations is granted if the following conditions are met for the parameters (θ, τ, σ):

σ > 0, τ > 0, θ > 0 (23a) γτ,σ ∈ [1; 2) (23b) γτ,σ > θ (23c) with γτ,σ := 2 -τ [1 -τ σλmax(B)] -1 L/2
, where L is the Lipschitz-continuity constant of ∇G, see Eq. [START_REF] Min | Fluorescent microscopy beyond diffraction limits using speckle illumination and joint support recovery[END_REF]. Within the convergence domain ensured by [START_REF] Denk | Two-photon laser scanning fluorescence microscopy[END_REF], the practical tuning of the parameter set (θ, τ, σ) is tedious as it may impair the convergence speed. We propose the following tuning strategy, which appeared to be very efficient. At first, we note that the step length τ relates only to the primal update (18a) whereas σ relates only to the dual update (18b) via ∆ (k) . In addition, the relaxation parameter θ scales both the primal and the dual steps [START_REF] Oh | Sub-Rayleigh imaging via speckle illumination[END_REF]. Considering only under-relaxation (i.e., θ < 1), (23c) is unnecessary and (23b) is equivalent to the following bound

σ ≤ σ with σ := 1/τ -L/2 λ -1 max (B). (24)
This relation defines an admissible domain for (τ, σ) under the condition θ < 1, see Fig. 10. Our strategy defines τ as the single tuning parameter of our PPDS iteration, the parameter σ being adjusted so that the dual step is maximized:

0 < τ < τ , σ = σ and θ = 0.99, (25)
with τ := 2/L. We set θ arbitrary close to 1 since practical evidence indicates that under-relaxing θ slows down the convergence rate. The numerical evaluation of the bounds τ and σ is application-dependent since they depend on L and λmax(B).

0 τ σ 1/L 2/L L 2 λ -1 max (B B)
Fig. 10. Admissible domain for (τ, σ) ensuring the global convergence of the PPDS iteration with θ ∈ (0; 1), see Equation [START_REF] Gu | Comparison of three-dimensional imaging properties between two-photon and single-photon fluorescence microscopy[END_REF].

C. Resolution of the joint Blind-SIM sub-problem

For our specific problem, the implementation of the PPDS iteration requires first the conjugate function [START_REF] Labouesse | Fluorescence blind structured illumination microscopy: a new reconstruction strategy[END_REF]: with h defined by (14b), the Fenchel conjugate is easily found and reads

P σh (ω) = vect (min {ωn, α}) . ( 26
)
The updating rule for the PPDS iteration then reads

q (k+1) ←-q (k) -θτ Bζ (k) (27a) ω (k+1) ←-ω (k) + θ∆ (k) (27b) with ∆ (k) = vect min{a (k) n , α} -ω (k) and a (k)
n the nth component of the vector a (k) defined in [START_REF] Golub | Matrix computation[END_REF]. We note that the positivity constraint is not enforced in the primal update (27a). Primal feasibility (i.e. positivity) therefore occurs only asymptotically thanks to the global convergence of the sequence [START_REF] Horstmeyer | Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM)[END_REF] toward the minimizer of the functional [START_REF] Rust | Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM)[END_REF]. Compared to FISTA, this behavior may be considered as a drawback of the PPDS iteration. However, we do believe that the ability of the PPDS iteration to "transfer" the hard constraint from the primal to the dual step is precisely the cornerstone of the acceleration provided by preconditioning. Obviously, such an acceleration requires that the preconditioner B is wisely chosen. For our joint Blind-SIM problem, the preconditioning matrix is derived from Geman and Yang semi-quadratic construction [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], [51, Eq. ( 6)]

B = 1 2 (C + β I d /a) -1 (28)
where I d is the identity matrix and a > 0 is a free parameter of the preconditioner. We choose C in the class of positive semidefinite matrix with a BCCB structure [START_REF] Vogel | Computational Methods for Inverse Problems[END_REF]Sec. 5.2.5]. This choice enforces that B is also BCCB, which considerably helps in reducing the computational burden: (i) B can be stored efficiently 12 and (ii) the matrix-vector product Bζ (k) in (27a) can be computed with O(N log N ) complexity by the bidimensional fast Fourier transform (FFT) algorithm. Obviously, if the observation model H is also a BCCB matrix built from the discretized OTF, the choice C = H t H in (28) leads to B = (∇ 2 g) -1 for a = 1. Such a preconditioner is expected to bring the fastest asymptotic convergence since it corrects the curvature anisotropies induced by the regular part g in the criterion [START_REF] Negash | Improving the axial and lateral resolution of three-dimensional fluorescence microscopy using random speckle illuminations[END_REF].

The PPDS pseudo-code for solving the joint Blind-SIM problem is given in Algorithm 1. This pseudo-code requires that L and λmax(B) are given for the tuning ( 25 // The primal step (Fourier domain)...

14 d (k) ←-ω (k) -2( h y -( γ + β) q (k) ) b; 15 q (k+1) ←-q (k) -θτ d (k) ; 16 // The dual step (direct domain)... 17 a (k) ←-FFT -1 ω (k) + σ( q (k) -2τ d (k) ) ; 18 ω (k+1) ←-(1 -θ) ω (k) + θ vect(min{a (k) n , α}); 19 // Prepare next PPDS iteration... 20 q (k) ←-q (k+1) ; ω (k) ←-FFT(ω (k+1) );
ρ ←-ρ + 1 M FFT -1 ( q (k) ) I0; 24 end
25 Final result: The joint Blind-SIM estimate is stored in ρ Algorithm 1: Pseudo-Code of the joint Blind-SIM PPDS algorithm, assuming that H is a BCCB matrix and C = H t H. The symbols and are the component-wise product and division, respectively. For the sake of simplicity, this pseudo-code implements a very simple stopping rule based on a maximum number of minimizing steps, see line 11. In practice, a more elaborated stopping rule could be used by monitoring the norm ||ζ (k) || defined by [START_REF] Jost | Optical sectioning and high resolution in single-slice structured illumination microscopy by thick slice blind-SIM reconstruction[END_REF] since it tends towards 0 as q (k) asymptotically reaches the constrained minimizer of the m-th nested problem. since H is rank deficient in our context, and the Lipschitz constant that reads L = λmax(B ∇ 2 g) can be further simplified as

L = a if a ≥ 1 ( γmax + β)( γmax + β/a) -1 otherwise, (29b)
with γmax the maximum of the square magnitude of the OTF components. From the pseudo-code, we also note that the computation of the primal update (27a) remains in the Fourier domain during the PPDS iteration, see line 14. With this strategy (possible because ∇g is a linear function), the computational burden per PPDS iteration 13is dominated by one single forward/inverse FFT pair, i.e., PPDS and FISTA have equivalent computational burden per iteration. We now illustrate the performance of the PPDS iterations for minimizing the penalized criterion involved in the joint Blind-SIM reconstruction problem shown in Fig. 9-(right). These simulations 9. The chosen initial-guess is q (0) = 0 for the primal variables and ω (0) = -∇g(q (0) ) for the dual variables. The preconditioning parameter is set to a = 1 and (θ, τ, σ) were set according to the tuning rule [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. For the sake of completeness, the curve of the FISTA iterations and the PDS iterations (i.e., the PPDS equipped with the identity preconditioning matrix B = I d ) are also reported.

f q (k) -f ( q)
were performed with a standard MATLAB implementation of the pseudo-code shown in Algorithm 1. We set a = 1 so that the preconditioner B is the inverse of the Hessian of g in [START_REF] Rust | Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM)[END_REF]. With this tuning, we expect that the PPDS iterations exhibit a very favorable convergence rate as long as the set of active constraints is correctly identified. Starting from initial guess q (0) = 0 (the dual variables being set accordingly to ω (0) = -∇g(q (0) ), see for instance [START_REF] Bertsekas | Nonlinear programming[END_REF]Sec. 3.3]), the criterion value of the PPDS iteration depicted in Fig. 11 exhibits an asymptotic convergence rate that can be considered as super-linear. Other tunings for a (not shown here) were tested and found to slow down the convergence speed. The pivotal role of the preconditioning in the convergence speed is also underlined since the PPDS algorithms becomes as slow as the standard proximal iteration when we set B = I d , see the "PDS" curve in Fig. 11. In addition, one can note from the reconstructions shown in Fig. 9 that the highfrequency components (i.e., the SR effect) are brought in the very early iterations. Actually, once PPDS is properly tuned, we always found that it offers a substantial acceleration with respect to the FISTA (or the standard proximal) iterates.

Finally, let us remind that the numerical simulations were performed with a BCCB convolution matrix H. In some cases, the implicit periodic boundary assumption 14 enforced by such matrices is not appropriate and a convolution model with a zero boundary assumption is preferable, which results in a matrix H with a BTTB structure. In such a case, the product of any vector by H t H can still be performed efficiently in O(N log N ) via the FFT algorithm, see for instance [START_REF] Vogel | Computational Methods for Inverse Problems[END_REF]Sec. 5.2.3]. This applies to the computation of ∇g(q (k) ) in the primal step [START_REF] Jost | Optical sectioning and high resolution in single-slice structured illumination microscopy by thick slice blind-SIM reconstruction[END_REF], according to [START_REF] Mukamel | Statistical deconvolution for superresolution fluorescence microscopy[END_REF]. In contrast, 14 Let us recall that the matrix-vector multiplication Hq with H a BCCB matrix corresponds to the circular convolution of q with the convolution kernel that defines H.

exact system solving as required by ( 21) cannot be implemented in O(N log N ) anymore if matrix H is only BTTB (and not BCCB). In such a situation, one can define C as a BCCB approximation of H t H, so that the preconditioning matrix B = (C +βI d ) -1 remains BCCB, while ensuring that B(H t H +βI d ) has a clustered spectrum around 1 as the size N increases [START_REF] Chan | Conjugate gradient methods for Toeplitz systems[END_REF]Th. 4.6].

Finally, another practical issue arises from the numerical evaluation of L. No direct extension of (29b) is available when H is BTTB but not BCCB. However, according to [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF], global convergence of the PPDS iterations is still granted if τ < 2/ L with L ≤ L. For instance, L := λmax(B) (||H||∞||H||1 + β) is an easy-to-compute upper bound of L.

V. CONCLUSION

The speckle-based fluorescence microscope proposed in [START_REF] Mudry | Structured illumination microscopy using unknown speckle patterns[END_REF] holds the promise of a super-resolved optical imager that is cheap and easy to use. The SR mechanism behind this strategy, that was not explained, is now properly linked with the sparsity of the illumination patterns. This readily relates joint Blind-SIM to localization microscopy techniques such as PALM [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF] where the image sparsity is indeed brought by the sample itself. This finding also suggests that "optimized" random patterns can be used to enhance SR, one example being the two-photon excitations proposed in this paper. Obviously, even with such excitations, the massively sparse activation process at work with PALM/STORM remains unparalleled and one may not expect a resolution with joint Blind-SIM greater than twice or three times the resolution of a wide-field microscope. We note, however, that this analysis of the SR mechanism is only valid when the sample and the illumination patterns are jointly retrieved. In other words, this article does not tell anything about the SR obtained from marginal estimation techniques that estimates the sample only, see for instance [START_REF] Min | Fluorescent microscopy beyond diffraction limits using speckle illumination and joint support recovery[END_REF]- [START_REF] Chaigne | Super-resolution photoacoustic fluctuation imaging with multiple speckle illumination[END_REF]. Indeed, the SR properties of such "marginal" techniques are rather distinct [START_REF] Idier | A theoretical analysis of the super-resolution capacity of imagers using unknown speckle illuminations[END_REF].

From a practical perspective, the joint Blind-SIM strategy should be tested shortly with experimental datasets. One expected difficulty arising in the processing of real data is the strong background level induced in the focal plane by the out-of-focus light. This phenomenon prevents the local extinction of the excitation intensity, hence destroying the expected SR in joint Blind-SIM. A natural approach would be to solve the reconstruction problem in its 3D structure, which is numerically challenging, but remains a mandatory step to achieve 3D speckle SIM reconstructions [START_REF] Negash | Improving the axial and lateral resolution of three-dimensional fluorescence microscopy using random speckle illuminations[END_REF]. The modeling of the out-of-focus background with a very smooth function is possible [START_REF] Orieux | Bayesian estimation for optimized structured illumination microscopy[END_REF] and will be considered for a fast 2D reconstruction of the sample in the focal plane.

Another important motivation of this work is the reduction of the computational time in joint Blind-SIM reconstructions. The reformulation of the original (large-scale) minimization problem is a first pivotal step as it leads to M sub-problems, all sharing the same structure, see Sec. II-A. The new preconditioned proximal iteration proposed in Sec. IV-B is also decisive as it efficiently tackles each sub-problem. In our opinion, this "preconditioned primal-dual splitting" (PPDS) technique is of general interest as it yields preconditioned proximal iterations that are easy to implement and provably convergent. For our specific problem, the criterion values are found to converge much faster with the PPDS iteration than with the standard proximal iterations (e.g., FISTA). We do believe, however, that PPDS deserves further investigations, both from the theoretical and the experimental viewpoints. This minimization strategy should be tested with other observation models and prior models. For example, as a natural extension of this work, we will consider shortly the Poisson distribution in the case of image acquisitions with low photon counting rates. The global and local convergence properties of PPDS should be explored extensively, in particular when the preconditioning matrix varies over the iterations. This issue is of importance if one aims at defining quasi-Newton proximal iterations with PPDS in a general context.

Fig. 1 .

 1 Fig. 1. [Row A] Lower-right quarter of the (160×160 pixels) groundtruth fluorescence pattern considered in [1] (left) and deconvolution of the corresponding wide-field image (right). The dashed (resp. solid) lines corresponds to the spatial frequencies transmitted by the OTF support (resp. twice the OTF support). [Row B] Positivity-constrained reconstruction from known illumination patterns: (left) M = 9 harmonic patterns and (right) M = 200 speckle patterns. The distance units along the horizontal and vertical axes are given in wavelength λ.

Fig. 2 .

 2 Fig. 2. [Row A] One product image qm = vect(ρn ×Im;n) built from one of the 200 illumination patterns used for generating the dataset: (left) a positive constant is added to the standard speckle patterns so that the lowest value is much greater that zero; (right) a positive constant is subtracted to the standard speckle patterns and negative values are set to zero. [Row B] Reconstruction of the product image qm that corresponds to the one shown above. [Row C] Final reconstruction ρ achieved with the whole set of illuminations -see Subsection II-B for details.

Fig. 3 .

 3 Fig. 3. Harmonic patterns: [Row A] One illumination pattern Im drawn from the set of regular (left) and distorted (right) harmonic patterns. [Row B] Corresponding penalized joint Blind-SIM reconstructions. [Row C] (left) Decreasing the number of phase shifts from 6 to 3 brings some reconstruction artifacts, see (B-left) for comparison. (right) Increasing the modulation frequency ||ν|| of the harmonic patterns above the OTF cutoff frequency prevents the super-resolution to occur. [Row D] Low-resolution image ym drawn from the dataset for a modulation frequency ||ν|| lying inside (left) and outside (right) the OTF domain-see Sec. III-A for details.

Fig. 4 .

 4 Fig. 4. Speckle patterns: [Row A] One speckle illumination such that NA ill = NA (left) and its "squared" counterpart (right). [Row B] Corresponding penalized joint Blind-SIM reconstructions from M = 1000 speckle (left) and "squared" speckle (right) patterns

Fig. 5 .

 5 Fig. 5. Speckle patterns (continued): Penalized joint Blind-SIM reconstructions from standard speckle (left) and "squared" speckle (right) patterns. The number of illumination patterns considered for reconstruction is M = 10 (A), M = 200 (B) and M = 10000 (C).

Fig. 6 .

 6 Fig. 6. Speckle patterns (continued): The correlation length of speckle and "squared" speckle patterns drives the level of super-resolution in the penalized joint Blind-SIM reconstruction: [Rows A] reconstruction from M = 10000 speckle patterns with NA ill = 0.5 NA (left) and from the corresponding "squared" random-patterns (right). [Rows B] idem with NA ill = 2 NA. [Rows C] idem with uncorrelated patterns.

Fig. 7 .

 7 Fig. 7. Processing of real and mock data: [Row A] Fluorescent beads with diameters of 100 nm are illuminated by 100 fully-developed (i.e., onephoton) speckle patterns through an illumination/collection objective (NA = 1.49). The sum of the acquisitions of the fluorescent light (left) and its Wiener deconvolution (middle) provide diffraction limited images of the beads. The joint Blind-SIM reconstruction performed with the hyper-parameters set to β = 5 × 10 -5 and α = 0.4 is significantly more resolved (right). The sampling rate used in these images is 32.5 nm, corresponding to an up-sampling factor of two with respect to the camera sampling. [Row B] STORM reconstruction of a marked rat neuron showing a lattice structure with a 190-nm periodicity (left). Deconvolution of the simulated wide-field image (middle). Joint Blind-SIM reconstruction of the sample obtained from 300 (one-photon) speckle patterns; the hyper-parameters are set to β = 2 × 10 -5 and α = 1.5 (right). The sampling rate of the STORM ground-truth image is 11.4 nm. The sampling rate of the joint Blind-SIM reconstruction is 28.5 nm, corresponding to an up-sampling factor of four with respect to the camera sampling. The distance units along the horizontal and vertical axes are given in wavelength λ coll , i.e., 520 nm in row A and 488 nm in row B.

Fig. 8 .

 8 Fig.8. Penalized Blind-SIM reconstructions from the dataset used to generate the super-resolved reconstruction shown in Fig.4(B-left). The hyper-parameter β was set to 10 -6 in any case, and α was set to 10 -3 (left) and 0.9 (right). For the sake of comparison, our tuning for the reconstruction shown in Fig.4(B-left) is β = 10 -6 and α = 0.3.

Fig. 9 .

 9 Fig. 9. Harmonic joint Blind-SIM reconstruction of the fluorescence pattern achieved by the minimization of the criterion (8) with 10, 50 or 1000 FISTA (abc) or PPDS (def) iterations. For all these simulations, the initial-guess is q (0) = 0 and the regularization parameters is set to (α = 0.3, β = 10 -6 ). The PPDS iteration implements the preconditioner given in[START_REF] Leterrier | Nanoscale architecture of the axon initial segment reveals an organized and robust scaffold[END_REF] with C = H t H and a = 1, see Sec IV-C for details.

13 6 ρ 11 /

 611 ): we get λmax(B) = 1/λmin(B -1 ) = a (2β) -1 (29a)12 Any BCCB matrix B reads H = F † ΛF with F the unitary discrete Fourier transform matrix, ' †' the transpose-conjugate operator, and Λ := Diag( b) where b := vect( bn) are the eigenvalues of B, see for instance [41, Sec. 5.2.5]. As a result, the storage requirement reduces to the storage of b. Given quantities:2 PSF h, Dataset {ym} M m=1 , Average intensity I0 ∈ R N + ; Regularization parameters: β, α ∈ R+; 4 PPDS parameters: a ∈ R+; θ ∈ (0, 1); τ ∈ (0, 2L); kmax ∈ N; ←-0; σ ←-σ [see[START_REF] Gu | Comparison of three-dimensional imaging properties between two-photon and single-photon fluorescence microscopy[END_REF]];7 h ←-FFT(h); γ ←-h * h; b ←-(2 γ + 2β/a);8 // The outer loop: processing each view ym... 9 for m = 1 • • • M do 10 y ←-FFT(ym); q (0) ←-FFT(q (0) m ); ω (0) ←-FFT(ω / The inner loop: PPDS minimization... 12 for k = 0 • • • kmax do 13

21 end 22 /

 2122 / Building-up the joint Blind-SIM estimate...

23

 23

Fig. 11 .

 11 Fig.11. Criterion value (upper plots) and distance to the minimizer (lower plots) as a function of the PPDS iterations for the reconstruction problem considered in Fig.9. The chosen initial-guess is q (0) = 0 for the primal variables and ω (0) = -∇g(q (0) ) for the dual variables. The preconditioning parameter is set to a = 1 and (θ, τ, σ) were set according to the tuning rule[START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. For the sake of completeness, the curve of the FISTA iterations and the PDS iterations (i.e., the PPDS equipped with the identity preconditioning matrix B = I d ) are also reported.

Whenever ρn = 0, the corresponding entry in the illumination pattern estimates (4b) can be set to Im;n = I 0;n /M for all m, hence preserving the positivity (2c) and the constraint (2b).

A constrained quadratic problem such as (3) is strictly convex if and only if the matrix H is full rank. In our case, however, H is rank deficient since its spectrum is the OTF that is strictly support-limited.

Assuming a fully-developed speckle, the fluctuation in Im;n is driven by an exponential pdf with parameter I 0 whereas the pdf of the "squared" pointwise intensity Jm;n := I 2 m,n is a Weibull distribution with shape parameter k = 0.5 and scale parameter λ = I 2 0 .

The MATLAB implementation of the PPDS pseudo-code Algorithm 1 requires less than 6 ms per iteration on a standard laptop (Intel Core M 1.3 GHz). For the sake of comparison, one FISTA iteration takes almost 5 ms on the same laptop.

ACKNOWLEDGMENTS

The authors are grateful to the anonymous reviewers for their valuable comments, and to Christophe Leterrier for the STORM image used in Section III.

Agence Nationale de la Recherche (ANR-12-BS03-0006