## Abstract

Sparsity constraint is a priori knowledge of the signal, which means that in some properly chosen basis only a small percentage of the total number of the signal components is nonzero. Sparsity constraint has been used in signal and image processing for a long time. Recent publications have shown that the Sparsity constraint can be used to achieve super-resolution of optical sparse objects beyond the diffraction limit. In this paper we present the quantum theory which establishes the quantum limits of super-resolution for the sparse objects. The key idea of our paper is to use the *discrete prolate spheroidal sequences* (DPSS) as the sensing basis. We demonstrate both analytically and numerically that this sensing basis gives superior performance of super-resolution over the Fourier basis conventionally used for sensing of sparse signals. The explanation of this phenomenon is in the fact that the DPSS are the eigenfunctions of the optical imaging system while the Fourier basis are not. We investigate the role of the quantum fluctuations of the light illuminating the object, in the performance of reconstruction algorithm. This analysis allows us to formulate the criteria for stable reconstruction of sparse objects with super-resolution. Our results imply that sparsity of the object is not the only parameter which describes super-resolution achievable via sparsity constraint.

© 2012 Optical Society of America

## 1. Introduction

Quantum Imaging [1, 2] is nowadays an established branch of quantum optics studying the ultimate performance limits of optical imaging allowed by the quantum nature of the light. One of the most important limits for any optical imaging scheme is that of its resolving capabilities. The classical diffraction limit established by Abbe and Rayleigh at the end of the 19th century, is nowadays recognized as a useful practical criterion but not a fundamental physical limit like, for example, the Heisenberg uncertainty relation. Numerous publications have demonstrated possibilities to improve the resolution beyond the Rayleigh limit. A natural question arises about the “real” quantum limit of optical super-resolution, i. e. resolution beyond the Rayleigh limit.

During the last ten years series of papers [3–6] have addressed the question of quantum limits of optical super-resolution in one- and two-dimensional coherent imaging for both continuous [3–5] and discrete [6] objects. In particular, the quantum theory developed in these papers, has established the standard quantum limit of super-resolution attained when the object is illuminated by a light wave in a multimode coherent quantum state. It was also shown that one can further increase the super-resolution factor beyond the standard quantum limit illuminating the object by the light in a multimode squeezed state.

In order to avoid confusion with other works on super-resolution, we would like to underline that by super-resolution we understand a possibility to exceed the Rayleigh diffraction limit imposed by a finite spatial bandwidth of the optical transfer function of the imaging system. Super-resolution in this sense is possible only when one has some a priori information about the object. The a priori information used in [3–6] is that the object has finite support. The spatial Fourier spectrum of such an object is an analytical function which allows for analytical continuation outside the transmission bandwidth of the optical imaging system. The standard quantum limit of super-resolution is set up by the quantum fluctuations of the optical field in the coherent state inside the support area of the object and the vacuum fluctuations outside this area.

In fact, the sparsity of the objects is another priori information that could be used for the super-resolution reconstruction. The mathematical theory about such signal processing was developed long time ago [7], and the practical applications of such method began to draw more attention since the concept of Compressed Sensing (CS) was proposed [8–11]. Recent publications have shown that the sparsity constraint can be used to achieve super-resolution in both conventional optical imaging [12, 13] and ghost imaging [14, 15]. However, to the best of our knowledge, the question of quantum limits of super-resolution achievable via sparsity constraint has not been addressed in the literature. In this paper we make a first step towards filling this gap. Namely, we combine the quantum theory of super-resolution from [3–6] with the additional a priori of sparsity in order to establish the quantum limits of super-resolution of sparse optical objects. Our analysis is an attempt to determine the ultimate performance capabilities via sparsity constraint imposed by the quantum fluctuations.

Let us mention that the problem of recovering sparse objects from band-limited measurements has a very long history and has been addressed in numerous publications. Compressed Sensing represents only one particular approach among many other approaches based on sparse and redundant representation theories and applications applied to image processing. We refer interested readers to a review article [7] and the references in this review.

The recovery algorithms mostly used in the processing of sparse signals come from the linear programming such as basis pursuit [10]. These algorithms use two different bases for representing the original object, known as the sensing and the representation basis. The first one is the basis in which the object is observed, while the second - in which the object is sparse. A typical example of such a pair is the Fourier basis and the canonical or spike basis. The key idea of our paper is to use *discrete prolate spheroidal sequences* (DPSS) as the sensing basis in the recovery algorithm via sparsity constraint. The DPSS were used recently for formulating the quantum theory of super-resolution of discrete optical objects [6]. Below we shall demonstrate that DPSS are an ideal tool for achieving super-resolution in reconstruction of discrete optical objects via sparsity constraint. Next, we investigate the role of the quantum fluctuations of the light illuminating sparse objects in the performance of the reconstruction algorithm. This analysis allows us to formulate the criteria of stable reconstruction of sparse objects with super-resolution. Our results imply that sparsity of the object is not the only parameter which characterizes super-resolution. Additional parameters specifying the distribution of the nonzero components of sparse objects are necessary in order to completely understand super-resolution via sparsity constraint.

The paper is organized as follows. In Sec. 2 we introduce the diffraction-limited optical imaging scheme of discrete optical objects and its eigenfunctions - discrete prolate spheroidal sequences, or DPSS. In Sec. 3 we describe the algorithm used for reconstruction of sparse objects via sparsity constraint. We demonstrate using numerical simulations that using DPSS as the sensing basis allows us to reconstruct noiseless sparse objects exactly, while the Fourier sensing basis produces errors in the reconstructed object even in the noiseless case. The quantum limits on super-resolution via the sparsity constraint are discussed in Sec. 4.

## 2. Optical imaging scheme and its eigenfunctions

In this section we shall describe the imaging scheme of discrete optical objects and its eigenfunctions - discrete prolate spheroidal sequences (DPSS), and discrete prolate spheroidal functions (DPSF). To simplify the description we shall start with classical analysis. We shall use the quantum theory in Sec. 4 where we address the question of quantum limits of super-resolution via the sparsity constraint.

In Fig. 1 we show a far-field optical imaging scheme of discrete objects. For simplicity, we consider a one-dimensional case, but our approach can be easily generalized to two-dimensional objects. A discrete object is formed by a sequence of rectangular pixels of size Δ within a finite support area *X*. The number *K* of pixels is given by *K* = *X*/Δ. As mentioned above, to achieve super-resolution one needs some a priori information about the object. In our case we know a priori that the object is confined within the area of size *X* and is identically zero outside. The spatial Fourier transform of such an object is an entire analytical function. Therefore, knowing the part of the spatial Fourier spectrum of the object within the diaphragm area in the Fourier plane, allows for analytical continuation of the spectrum and thus, super-resolution.

A thin lens *L* placed at the focal distance from the object pane, performs the spatial Fourier transform of the object onto the Fourier plane. The diaphragm of size *d* in the Fourier plane introduces the finite spatial transmission bandwidth of the imaging system (we neglect diffraction on the lens *L*). Only the spatial Fourier components of the object within the aperture are transmitted through the system and can be detected in the Fourier plane. We shall assume that these transmitted Fourier components are detected by means of a homodyne detection with a local oscillator and a sensitive CCD camera, so that both quadrature components of the field can be measured.

Essentially, the scheme in Fig. 1 is the same imaging system as that considered in Ref. [6]. The only notable difference is that here, contrary to Ref. [6], we are dealing with sparse discrete objects, which means that only a part of the total number *K* of pixels has nonzero values. In an example shown in Fig. 1 we have only two nonzero pixels out of total five. In general, we shall assume that the number of nonzero pixels is equal to *P*. The ratio *P/K* is called sparsity parameter in the CS theory [8–11]. Even if the positions of these *P* nonzero pixels are unknown as well as the values of the corresponding field amplitudes, sparsity gives additional a priori information which should further increase the super-resolution factor over the conventional dense discrete objects studied in Ref. [6].

Let us introduce the dimensionless spatial coordinates in the object plane as *s* = 2*x*/*X*, and in the Fourier plane as *ξ* = 2*y/d*. In these dimensionless coordinates both the support area of the object and the area of the aperture are normalized to [−1, 1]. The complex field amplitudes *a*(*s*) in the object plane, and *f*(*ξ*) in the Fourier plane are related by the Fourier transform performed by the lens *L*, which in terms of these dimensionless variables reads,

In the dimensional coordinates *s* the size of a pixel becomes Δ = 2/*K*. For symmetry reasons we will consider the case of odd number *K* = 2*M* +1. Let us introduce the “pixellized” complex amplitudes *a*(*m*),

*m*= 0, ±1, ±2,... is the position of the center of

*m*-th pixel. The support area corresponds to

*m*= 0, ±1,... ±

*M*.

As explained in Ref. [6], one can introduce a coarse-grained complex amplitude of the object constructed from the “pixelized” complex amplitudes *a*(*m*) and the pixel frame function. The spatial Fourier spectrum *F*(*ξ*) of such a coarse-grained complex amplitude is a product of two terms: *f*(*ξ*) containing the information about the coefficients *a*(*m*), and *P*(*ξ*) which is the Fourier spectrum of the pixel frame function. Since the latter term does not carry information about the object coefficients, we shall assume that this part can be eliminated from the measured spectrum *F*(*ξ*). Therefore, in what follows we shall deal with the corrected Fourier spectrum *f*(*ξ*) = *F*(*ξ*)/*P*(*ξ*) given by,

The discrete prolate spheroidal sequences (DPSS) *ψ _{μ}* (

*m*),

*μ*= 0, 1...2

*M*, are the

*K*-dimensional eigenvectors satisfying the following equation [16]:

*λ*are the nondegenerate real eigenvalues such that 1 >

_{μ}*λ*

_{0}> ... >

*λ*

_{2M}> 0. The DPSS have the following properties:

*ξ*| < 1/2

*W*. Let us note that our dimensionless Fourier frequencies

*ξ*differ from the Fourier frequencies mostly used in the literature on the DPSS by a factor 1/

*W*.

Since the Fourier spectrum given by Eq. (3) is a continuous function, we need another basis set in the Fourier plane. This set is given by the discrete prolate spheroidal functions (DPSF) Ψ* _{μ}* (

*ξ*), the eigenfunctions of the following equation:

*λ*,

_{μ}*μ*= 0, 1...2

*M*as in Eq. (4). The DPSF are doubly orthogonal:

*μ*= 0, 1...2

*M*.

Using DPSS *ψ _{μ}* (

*m*),

*μ*= 0, 1...2

*M*we can introduce the complete orthonormal basis sequences in the support area of the object as,

*a*(

*m*) in the object plane, and the corresponding Fourier coefficients

*f*(

*ξ*) in the Fourier plane, As demonstrated in Ref. [6] the classical complex coefficients

*a*in the object plane and

_{μ}*f*in the Fourier plane are related by,

_{μ}Let us assume that we can measure the classical Fourier amplitudes *f*(*ξ*) within the transmission area of the diaphragm using a homodyne detection scheme and a sensitive CCD camera with the size of pixels chosen appropriately in order to resolve the features in the continuous Fourier spectrum. Then we can calculate the coefficients
${a}_{\mu}^{(r)}$ of the reconstructed object as,

## 3. Super-resolution of sparse noiseless objects via sparsity constraint

Most natural signals have only a small number of nonzero components when they are expressed in an appropriately chosen basis. An example is a digital image taken by a megapixel digital camera. Although typically almost all the image pixels have nonzero values, most of digital images turn out to be sparse when expressed in an appropriately chosen wavelet basis. Such signals possess a priori information called “sparsity”. This sparsity property of digital images explains why one can effectively compress most of them almost without loosing their information content.

In the reconstruction processing via sparsity constraint, there are two different basis sets for representing an object: a sensing basis *ϕ _{μ}* (

*m*) and a representation basis

*ϕ*′

*(*

_{μ}*m*). The first one is the basis in which the object is observed, while the second one is the basis in which the object is sparse. We shall assume that these two basis sets are orthonormal, however this assumption is not necessary in general. Applying this terminology to our discrete optical object, we shall write its complex amplitudes

*a*(

*m*) in the representation basis as,

*c*′

*and the sensing coefficients*

_{μ}*c*. For sparse objects only a part of representation coefficients

_{μ}*c*′

*, usually small, is nonzero. In this paper we shall assume that the object is sparse in the spike basis*

_{μ}*ϕ*′

*(*

_{μ}*m*) =

*δ*

_{μ,m+M+1},

*μ*= 1, …,

*K*= 2

*M*+ 1;

*m*= −

*M*,…,

*M*, and that only

*P*≪

*K*representation coefficients

*c*′

*are nonzero.*

_{μ}Using Eqs. (15) and (16) it is straightforward to obtain the following linear relationship between the sensing and the representation coefficients,

where the matrix*A*with components

*A*, is known in the theory of compressed sensing as the sensing matrix.

_{μν}If we do not have a priori information about the sparsity of the signal *a*(*m*), all the *K* representation coefficients *c*′* _{ν}*,

*ν*= 1,...,

*K*, are nonzero. Then, in order to solve the linear system Eq. (17), we need to know all the

*K*sensing coefficients

*c*,

_{ν}*ν*= 1,...,

*K*. Thus we obtain a square

*K*×

*K*sensing matrix

*A*. The situation is drastically different for sparse signals when we know a priori that only

*P*≪

*K*representation coefficients

*c*′

*are nonzero, even if we do not know which ones. By using the priori information of sparsity to solve Eq. (17), we can use only*

_{ν}*Q*<

*K*sensing coefficients. This statement seems paradoxical, since the sensing matrix now becomes rectangular

*Q*×

*K*, and Eq. (17) seems to be underdetermined, because we have only

*Q*<

*K*equations for

*K*unknowns

*c*′

*. Nevertheless, the priori information of sparsity arrives at recovering the signal using an additional constraint of*

_{ν}*l*

_{1}-norm minimization [8–11] for the solution of Eq. (17).

Another important notion in addition to the notion of sparsity is that of the coherence between the sensing basis Φ = {*ϕ _{μ}*} and representation basis Φ′ = {

*ϕ*′

*}. The coherence*

_{μ}*γ*between these two bases is defined as [10]:

*ϕ*,

_{μ}*ϕ*′

*〉 the scalar product,*

_{ν}*γ*is directly related to the number of measurements

*Q*required for the reconstruction of a sparse object as

*Q*∼

*γ*

^{2}

*P*log

*K*. From this relationship it becomes clear why reconstruction via sparsity constraint deals mainly with low-coherence basis pairs Φ and Φ′. An example of such a pair is the canonical or spike basis

*ϕ*′

*(*

_{μ}*m*) =

*δ*

_{μ,m+M+1}and the Fourier basis ${\varphi}_{\mu}(m)=\text{exp}\left[-2\pi i\mu m/K\right]/\sqrt{K}$, which realize the maximum incoherence,

*γ*= 1. This is the reason why we have chosen the Fourier basis for comparison with our results obtained using DPSS as the sensing basis.

Before applying the sparsity constraint for reconstruction of an original sparse object using DPSS as the sensing basis, we have evaluated numerically the coherence *γ* between the representation basis Φ′, which is the spike basis *ϕ*′* _{μ}*(

*m*) =

*δ*

_{μ,m+M+1},

*μ*= 1,...,

*K*= 2

*M*+ 1;

*m*= −

*M*,...,

*M*in our model, and the sensing basis Φ, that we choose as the orthonormal DPSS

*φ*(

_{μ}*m*) defined in Eq. (9):

*ϕ*(

_{μ}*m*) =

*φ*

_{μ}_{+1}(

*m*). For completeness we have also evaluated another coherence for the case when the representation basis Φ′ is the Fourier basis, i. e. ${\varphi}_{\mu}^{\prime}(m)=\text{exp}\left[-2\pi i\mu m/K\right]/\sqrt{K}$. To distinguish these two cases we have denoted the latter coherence as Γ. This coherence may be useful for reconstruction of objects which are

*sparse in the Fourier representation basis*. However, we will not consider this case in the present paper. The results of our numerical evaluation of both coherences,

*γ*and Γ, are shown in Fig. 2. We have plotted

*γ*and Γ as functions of the pixel number

*K*and the parameter 0 ≤

*W*≤ 1/2 which determines the spatial bandwidth of

*φ*(

_{μ}*m*) in the dimensional coordinates

*ξ*[see Eq. (7)].

In Fig. 2 dots indicate the grid at which the values of the coherences were numerically evaluated, while the continuous lines are the fitting curves. As follows from (a) and (c), both coherences *γ* and Γ for fixed *W* are growing linearly with
$\sqrt{K}$. The speed of this linear growth remains smaller than unity for any value of *W*. One can notice that the coherence *γ* between the spike and the DPSS basis is always lower than the coherence Γ between the Fourier and the DPSS basis. From (b) and (d) one can observe the symmetry of *γ* and Γ as functions of *W* with respect to *W* = 1/4. The reason for this symmetry comes from the symmetry property of DPSS [16]:

We turn now to presentation of our results for reconstruction with super-resolution of noiseless sparse objects *a*(*m*) via sparsity constraint. Following Ref. [10], we assume that we can measure sensing coefficients *c _{μ}* in the sensing basis

*ϕ*(

_{μ}*m*) with

*μ*= 1, 2,...,

*Q*<

*K*. We shall use two different sensing bases, the Fourier basis and the DPSS basis. Here it is important to stress the difference between the sensing assumption in the standard CS theory and in our case, using the imaging scheme from Fig. 1. In the CS theory it is assumed that we have access to the whole spectrum

*f*(

*ξ*), for instance, using the Fourier sensing basis we obtain the sensing coefficients

*c*by sampling the whole spectrum randomly or regularly at

_{μ}*Q*discrete points. In our case we do not have access to the whole spectrum, but only to its low-frequency part within the diaphragm |

*ξ*| ≤ 1. Therefore, we shall assume that using the Fourier sensing basis we obtain the sensing coefficients

*c*by sampling the Fourier spectrum

_{μ}*f*(

*ξ*) at equidistant points

*ξ*,

_{μ}*μ*= 1, 2,...,

*Q*within the diaphragm |

*ξ*|

*≤*1:

*c*=

_{μ}*f*(

*ξ*).

_{μ}Using DPSS as the sensing basis corresponds to projection of the spatial Fourier spectrum *f*(*ξ*) onto the continuous DPSF Φ* _{μ}* (

*ξ*) defined in Eq. (10) within the diaphragm area,

*c*are obtained as ${c}_{\mu}={f}_{\mu}/\sqrt{{\lambda}_{\mu}}$. We firstly arrange the DPSF and DPSS such that the corresponding eigenvalues are in a descending order 1 >

_{μ}*λ*

_{0}> ... >

*λ*

_{2M}> 0, and then gradually increase the number of the DPSS and DPSF used for reconstruction in such arranged order until the reconstructions succeed. Let us note that practically the measurement corresponding to Eq. (22) would require discretization of the continuous integral in this equation and using of a sensitive CCD camera with the size of pixels chosen appropriately in order to resolve the features in the continuous Fourier spectrum. This type of measurement can be thought of as wasteful in the spirit of CS. Therefore, we would like to stress that our goal here is not minimization of the number of measurements as in CS, but achieving super-resolution in the reconstructed image via sparsity constraint.

Having the sensing coefficients *c _{μ}* we attempt to reconstruct the original object

*a*(

*m*) using the

*l*

_{1}-minimization in the spirit of sparsity constraint. Precisely, we are looking for a solution of the following convex optimization program,

*A*is the sensing matrix corresponding to each sensing basis, Fourier and DPSS. In order to solve Eq. (23) we employ the basis pursuit method [18].

The results of reconstruction of a sparse object *a*(*m*) with super-resolution via sparsity constraint are presented in Fig. 3. We choose a discrete sparse object *a*(*m*) with *P* = 5 nonzero pixels out of total *K* = 51. The positions of these 5 pixels are selected randomly and their amplitudes can take positive or negative value within the interval [−1, 1]. In three columns (1)–(3) we take three different values of the space-bandwidth product *c*, equal to 30, 20, and 10. This parameter determines the relative size of the diaphragm with respect to the period of the spatial Fourier spectrum of the object, as shown in the row (a). In the row (b) we plot by the black dots the original sparse object *a*(*m*), and by the red circles the corresponding diffraction-limited images that one would observe in the image plane of the scheme from Fig. 1. In the row (c) we show the original object *a*(*m*) (black dots) together with reconstructed object (red circles) when we use the Fourier basis as the sensing basis. And in the last row (d) we present the original object together with the reconstructed object using the DPSS basis as the sensing basis.

Comparing the two last rows (c) and (d) one can immediately see that the DPSS sensing basis provides a superior performance in the reconstruction of the original object *a*(*m*) over the Fourier sensing basis. Indeed, while the reconstructed object in the row (c) contains errors, the reconstruction with the DPSS sensing basis reproduces the original object exactly for all three values of the space-bandwidth product *c*. In order to characterize quantitatively the reconstruction errors for the Fourier sensing basis we have evaluated the following relative error,

*a*

^{(r)}(

*m*) is the reconstructed object obtained using the Fourier sensing coefficients ${c}_{\mu}={f}_{\mu}/\sqrt{{\lambda}_{\mu}}$. For three values of the space-bandwidth product

*c*,

*c*= 30, 20, and 10, in three columns (1)–(3) of Fig. 3, we have obtained the relative errors of reconstruction 0.091, 0.845, and 0.858, respectively.

It might seem counterintuitive that the DPSS sensing basis provides better performance than the Fourier sensing basis. Indeed, the spike-DPSS basis pair has higher coherence than the spike-Fourier pair and should, therefore, correspond to higher number of measurements for reconstruction of the original object according to the CS theory. Let us recall, however, that the problem of super-resolution is different from traditional reconstruction of sparse objects via CS method from a limited number of measurements. The difference is that in the CS scenario one has an access to the whole Fourier spectrum, while in the super-resolution scenario we have only a low-pass part of the Fourier spectrum available. Therefore, the conclusions of the CS theory about the total number of measurements are not applicable in this case. This is illustrated in Fig. 3 where the spike-Fourier pair does not provide exact reconstruction while the spike-DPSS does at expense of increasing the number of basis functions.

Let us explain now why we claim achieving super-resolution in the reconstructed object *a*^{(r)} (*m*) in the last row of Fig. 3. As explained in Ref. [6] the meaning of super-resolution for discrete objects (including sparse discrete objects) is different from continuous objects. Indeed, in the discrete case an object is completely characterized by its *K* pixel components and, therefore, super-resolution cannot be thought of as an attempt to resolve the fine details of the object, as in the continuous case. However, both discrete and continuous situations are closely related in terms of the spatial Fourier spectra. The spatial Fourier spectrum of a discrete object with *K* components is a periodic function with a period |*ξ*| = 1/2*W* in terms of the dimensionless Fourier frequency *ξ*. Since the size of the diaphragm in Fig. 1 is given by |*ξ*| = 1, for *W* = 1/2 all the Fourier components pass through the diaphragm and the object can be perfectly reconstructed from its Fourier components within the diaphragm. There is no super-resolution in this case. For *W* < 1/2 part of spectral components of the object is absorbed by the diaphragm. Super-resolution is an attempt to recover these absorbed components using the analyticity of the spectrum of an object with finite support. As demonstrated in Ref. [6], the super-resolution factor (SF) of the scheme in Fig. 1 is given by the ratio of the diagraph size to the spectrum period: *SF* = 1/2*W*. For *K* = 51 and *c* = 30, 20, 10, we have *W* = *c*/(*πK*) = 0.187, 0.125, 0.062 and the super-resolution factor *SF* = 2.7, 4, 7, respectively.

When we use the Fourier sensing basis, to obtain reconstruction through the basis pursuit method for *c* = 30, we have found empirically that we need at least *Q* = 23 equidistant sampling points within the diaphragm in the Fourier plane in Fig. 1. This number roughly corresponds to the number of linearly independent discrete Fourier components *Q*̃ inside the aperture, given by *Q*̃ = 2*WK*. For *c* = 30 we have *Q*̃ = 19. We have observed that as long as the number of measurements is larger than the number of linearly independent discrete Fourier components within the diaphragm, *Q* ≥ *Q*̃, increasing the number of measurements *Q* used in the basis pursuit method does not improve the quality of the reconstructed object. This can be explained by the fact that increasing the number of measurements *Q* over the limit *Q* = *Q*̃ corresponds to sub-sampling within the diaphragm. Such sub-sampling gives additional Fourier components which are linearly dependent with *Q*̃ Fourier components, but it does not help in recovering *K* − *Q*̃ Fourier components which were absorbed. This explanation is confirmed by the reconstructed objects in columns (2) and (3), row (c), where we use the same number of measurements *Q* = 23 for reconstruction in order to demonstrate deterioration of the reconstructed object with decreasing space-bandwidth product c. Although the number of measurements *Q* here is larger than the number of linearly independent discrete Fourier components within the diaphragm *Q*̃: *Q*̃ = 13 for *c* = 20, and *Q*̃ = 8 for *c* = 10, it is not possible to obtain exact reconstructions.

Reconstruction with DPSS sensing basis is different from that using the Fourier sensing basis. Now the sensing coefficients *c _{μ}* are obtained from projections of the Fourier spectrum

*f*(

*ξ*) within the diaphragm onto the continuous DPSF following Eq. (22). Using the theory from Ref. [6], one can show that in this case the sensing coefficients obtained using only part of the spectrum within the diaphragm, and those obtained using the total spectrum, are linearly dependent. That is why increasing the number Q of DPSS basis for the reconstruction we can always obtain perfect reconstruction of the original object in the noiseless case.

In the examples shown in the last row (d) for *c* = 30, 20, 10, we used the first Q = 27,28,30 DPSF and DPSS with the largest eigenvalues, respectively. These are the lowest numbers of sensing coefficients which guarantee stable reconstruction through the basis pursuit algorithm. In conclusion of this Section let us note that our results are not restricted to the particular object shown in Fig. 3. Indeed, we have repeatedly obtained perfect reconstruction of different random objects with DPSS sensing basis in the noiseless case, while there is no such guarantee with the Fourier sensing basis.

## 4. The role of the quantum fluctuations in super-resolution via sparsity constraint

In previous section we have shown that using DPSS as the sensing basis a sparse noiseless object can be exactly reconstructed through convex optimization to achieve super-resolution. In this section we shall address the question of the role of the quantum fluctuations in the reconstruction via sparsity constraint. As we shall demonstrate below, quantum fluctuations impose the ultimate performance limit on super-resolution in recovering discrete sparse objects similarly to the case of super-resolution of discrete dense objects considered in Ref. [6]. Since the quantum theory for the imaging scheme in Fig. 1, developed in Ref. [6], is valid for sparse objects as well, we shall give here its brief overview.

In the quantum theory the classical complex amplitudes *a*(*m*) and *a*^{★}(*m*) become the photon annihilation and creation operators *â*(*m*) and *â*^{†}(*m*) satisfying the canonical commutation relations,

*f*̂(

*ξ*),

*f*̂

^{†}(

*ξ*) in the Fourier plane obey the commutation relations in the continuous coordinate

*ξ*,

*L*in Fig. 1, remains valid for the quantum-mechanical operators, Using DPSS

*ψ*(

_{μ}*m*) as basis functions in the object plane and DPSF Ψ

*(*

_{μ}*ξ*) in the Fourier plane, we can write the Fourier transform Eq. (27) as a unitary transformation of the photon annihilation operators. To obtain this unitary transformation we split the object plane into the “core” region |

*m*| ≤

*M*, corresponding to the object support area, and the “wings” region |

*m*| >

*M*, outside the support area. The orthonormal basis sequences in these two regions are given by,

*ξ*| ≤ 1 corresponding to the transmission area of the diaphragm, and the wings region |

*ξ*| > 1 of the absorption area. The orthonormal basis functions in these regions are,

*â*(

*m*) in the object plane as,

*f*̂(

*ξ*) in the Fourier plane,

*â*and

_{μ}*f*̂

*are the annihilation operators of the discrete prolate modes in the corresponding core regions, and*

_{μ}*b*̂

*and*

_{μ}*ĝ*- the annihilation operators on the wings regions. The terms

_{μ}*â*

_{⊥}(

*m*) and

*f*̂

_{⊥}(

*ξ*) are added in order to satisfy the commutation relations Eq. (25) and Eq. (26). These terms are orthogonal to the prolate modes and will not appear in the unitary transformation.

Substituting Eqs. (30) and (31) into Fourier transform Eq. (27) and using Eqs. (8), (28)–(29) we arrive at the following unitary transformation:

*b*̂

*. The operator-valued coefficients ${\widehat{a}}_{\mu}^{(r)}$ of the reconstructed object become,*

_{μ}*b*̂

*describing the vacuum fluctuations in the wings area. These vacuum fluctuations do not allow us to reconstruct the high-order modes*

_{μ}*â*, because the eigenvalues

_{μ}*λ*become rapidly very small after the index

_{μ}*μ*has attained its critical value (see below). The case of coherent input state in the core region and the vacuum state in the wings region determines the

*standard quantum limit of super-resolution*for discrete dense objects [6].

For numerical simulations of the quantum fluctuations we choose a *c*-number representation of the quantum-mechanical operators *â _{μ}* and

*b*̂

*in Eq. (33) corresponding to the*

_{μ}*antinormal*ordering [19]. In this representation the operators

*â*and

_{μ}*b*̂

*become the*

_{μ}*c*-number Gaussian stochastic variables

*α*and

_{μ}*β*,

_{μ}*a*= 〈

_{μ}*â*〉 is the mean value of the annihilation operator in the core region. It coincides with the classical decomposition coefficients in Eq. (11). The mean values of the annihilation operators in the wings area are zero, 〈

_{μ}*b*̂

*〉 = 0, corresponding to zero classical amplitudes outside the support area of the object. The quantum fluctuations are accounted for by stochastic Gaussian fluctuations*

_{μ}*δα*and

_{μ}*δβ*. We introduce the quadrature components of these fluctuations as follows:

_{μ}*m*| ≤

*M*and the vacuum state in the wings area |

*m*| >

*M*of the object plane the nonzero correlations of the quadrature fluctuations are given by,

In order to understand the role of the quantum fluctuations in the reconstruction via sparsity constraint, we have to recall the behavior of the eigenvalues *λ _{μ}* of DPSS with index

*μ*. As already mentioned above, all the eigenvalues

*λ*are real and such that 1 >

_{μ}*λ*

_{0}> ... >

*λ*

_{K−1}> 0. An important parameter, characterizing the behavior of the eigenvalues

*λ*is the Shannon number

_{μ}*S*, defined as

*S*= 2

*c*/

*π*. For

*μ*≤

*S*the eigenvalues

*λ*are close to unity, while for

_{μ}*μ*>

*S*they rapidly approach zero with growing

*μ*[16, 20].

Using these properties of the eigenvalues *λ _{μ}* and Eq. (33) for the reconstructed operator-valued coefficient
${\widehat{a}}_{\mu}^{(r)}$ we can qualitatively understand the role of the quantum fluctuations in reconstruction of the original object. As follows from Eq. (33), to obtain the reconstruction which is stable with respect to the vacuum fluctuations, we have to guarantee that the number

*Q*of DPSS used as the sensing basis, is not larger than the Shannon number of the optical system,

*Q*≤

*S*.

To verify this behavior, we have performed numerical simulations of the reconstruction of a sparse object taking into account the quantum fluctuations, for two different cases, *Q* ≤ *S* and *Q* > *S*. We have considered the case of the space-bandwidth product *c* = 20 which corresponds to the Shannon number *S* = 12.7. We have used two different sparse objects with the same total number of pixels, *K* = 51, and the number of nonzero pixels *P* = 5. For the first object the number *Q* of sensing coefficients necessary for achieving reconstruction via the basis pursuit method, is *Q* = 13, while for the second object this number is *Q* = 19. It should be noted that the number *Q* of sensing coefficients needed for the reconstruction of a sparse object depends on the particular distribution of *P* nonzero pixels over the total number of pixels *K*.

We would like to underline here that the number of sensing coefficients needed for successful reconstruction of an object via sparsity constraint depends not only on the sparsity of the object, but also on the particular distribution of nonzero pixels. To the best of our knowledge, at present there is no quantitative relationship between the number of measurements and the positions of the nonzero pixels. Therefore we have tried several different sparse objects before arriving at these two particular examples in Fig. 4. There is a general belief that “when the positions of *K* nonzero pixels are widely scattered or at least with the same sign of phase, the reconstruction via sparsity constraint is more likely to succeed” [12]. We refer interested readers to Refs. [21] and [22] for more details.

For the numerical simulations of the quantum fluctuations, we follow the *c*-number representation of the operators *a _{μ}* and

*b*in Eq. (34), where the quantum fluctuations are accounted for by stochastic Gaussian fluctuations

_{μ}*δα*and

_{μ}*δβ*as the combination of two quadrature components in Eq. (35). The relative role of the quantum fluctuations in the reconstruction process depends on the signal-to-noise ratio (SNR) of the scheme,

_{μ}*N*̂〉 is the mean number of photons passed through the object area during the observation time [4], and 〈(Δ

*N*̂)

^{2}〉 is its variance. For an object in a coherent state and the vacuum state outside the object this SNR is equal to the mean number of photons 〈

*N*̂〉,

In Fig. 4 we present the results of our simulations. The rows (a)–(c) correspond to the case *Q* = 13 while the rows (d)–(f) to the case *Q* = 19. The values of the SNR are: 10^{12} in (a) and (d), 10^{8} in (b) and (e), and 10^{4} in (c) and (f). For each value of the SNR we present 10 random realizations of the reconstructed object, indicated by the circles of different color. Figure 4 confirms our criterion of stable reconstruction according to the number *Q* of sensing coefficients with respect to the Shannon number *S*. Indeed, for the first object with *Q* ≤ *S* we obtain stable reconstructions in large range of SNRs, while for the second object with *Q* > *S* the reconstruction becomes unstable for lower values of the SNR.

As we have already mentioned, the number *Q* of sensing coefficients needed for stable reconstruction depends on a particular distribution of *P* nonzero pixels over the total number of *K* pixels. We have performed statistical analysis of stable reconstructions of different objects with the same sparsity parameter *P*/*K* but different patterns of nonzero components. For each fixed ratio *P*/*K* we have chosen 100 different objects with random positions of nonzero pixels and random amplitudes *a*(*m*) from the interval between −1 and 1. We have fixed the number *Q* of sensing coefficients by requiring the smallest eigenvalue *λ _{Q}* in the decomposition over DPSS to be

*λ*= 0.01.

_{Q}The results of our statistical analysis are summarized in Table 1 which shows the percentage of stably reconstructed objects for different values of *P*, *K*, and different super-resolution factors *SF*. As follows from this table, one can obtain higher super-resolution factors for smaller values of the sparsity parameter *P*/*K*. This result corroborates the key role of a priori information in super-resolution. In this paper we do not study in detail the super-resolution factor as function of the sparsity parameter *P*/*K* and the SNR. This study will be the subject of a separate publication elsewhere.

In conclusion, we have presented the quantum theory of super-resolution of optical sparse objects via sparsity constraint. Our analysis demonstrates that using discrete prolate spheroidal sequences or DPSS as sensing basis in the reconstruction algorithm, one can obtain exact reconstruction of original noiseless sparse objects with high super-resolution factor. In our numerical simulations we have observed better performance of the DPSS-spike basis pair over the Fourier-spike pair. The latter basis pair does not allow for exact reconstruction of noiseless objects with high super-resolution factors. Next, we have applied the quantum theory of super-resolution to the scenario of the reconstruction via sparsity constraint. We have formulated a qualitative criterion for stable reconstruction of sparse objects in terms of the number *Q* of sensing coefficients and the Shannon number *S* of the optical imaging system. In order to achieve stable reconstruction with respect to both quantum and classical fluctuations one needs to guarantee that *Q* ≤ *S*. Our numerical simulations corroborate this criterion.

## Acknowledgments

We would like to thank Carlos Fernandez-Granda for his help with numerical simulations. Mikhail I. Kolobov is Fulbright Visiting Scholar at the Department of Mathematics and Statistics, Stanford University. He expresses his gratitude to Emmanuel Candès for hospitality and stimulating discussions during his stay at Stanford. He also acknowledges financial support of the Fulbright Commission and the French National Agency (ANR) through Nanoscience and Nanotechnology Program (Project NATIF n°ANR-09-NANO-012-01).

## References and links

**1. **M. I. Kolobov (Ed.), *Quantum Imaging* (Springer, NY, 2006).

**2. **S. J. Bentley, *Principles of Quantum Imaging: Ghost Imaging, Ghost Diffraction, and Quantum Lithography* (Taylor & Francis, Boca Raton, 2010).

**3. **M. I. Kolobov and C. Fabre, “Quantum limits on optical resolution,” Phys. Rev. Lett. **85**, 3789–3792 (2000). [CrossRef] [PubMed]

**4. **V. N. Beskrovnyy and M. I. Kolobov, “Quantum limits of super-resolution in reconstruction of optical objects,” Phys. Rev. A **71**, 043802 (2005). [CrossRef]

**5. **V. N. Beskrovnyy and M. I. Kolobov, “Quantum-statistical analysis of superresolution for optical systems with circular symmetry,” Phys. Rev. A **78**, 043824 (2008). [CrossRef]

**6. **M. I. Kolobov, “Quantum limits of superresolution for imaging discrete subwavelength structures,” Opt. Express **16**, 58–66 (2008). [CrossRef] [PubMed]

**7. **M. Elad, M. A. T. Figueiredo, and Y. Ma, “On the Role of Sparse and Redundant Representations in Image Processing,” Proc IEEE **98**, 972–982 (2010). [CrossRef]

**8. **E. J. Candès, J. Romberg, and T. Tao,“Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory **52**, 489–509 (2006). [CrossRef]

**9. **E. J. Candès and T. Tao, “Near-optimal signal recovery from random projections: Universal encoding strategies?” IEEE Trans. Inf. Theory **52**, 5406–5425 (2006). [CrossRef]

**10. **E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. **25**, 21–30 (2008). [CrossRef]

**11. **D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory **52**, 1289–1306 (2006). [CrossRef]

**12. **S. Gazit, A. Szameit, Y. C. Eldar, and M. Segev, “Super-resolution and reconstruction of sparse sub-wavelength images,” Opt. Express **17**, 23920–23946 (2009). [CrossRef]

**13. **Y. Shechtman, S. Gazit, A. Szameit, Y. C. Eldar, and M. Segev, “Super-resolution and reconstruction of sparse images carried by incoherent light,” Opt. Lett. **35**, 1148–1150 (2010). [CrossRef] [PubMed]

**14. **W. Gong and S. Han, “Super-resolution ghost imaging via compressive sampling reconstruction,” arXiv [0910.4823v1] (2009).

**15. **W. Gong and S. Han, “Super-resolution and reconstruction of far-field ghost imaging via sparcity constraints,” in the *Proceedings of SPARS’11-Singal Processing with Adaptive Sparse Structured Representations*, (Edinburgh, Scotland, 2011), p. 91. [PubMed]

**16. **D. Slepian, “Prolate spheroidal wave functions, Fourier analysis and uncertainty-V: The discrete case,” Bell System Tech. J. **57**, 1371–1430 (1978).

**17. **D. L. Donoho and X. Huo, “Uncertainty principles and ideal atomic decomposition,” IEEE Trans. Inf. Theory **47**, 2845–2862 (2001). [CrossRef]

**18. **S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Review **43**, 129–159 (2001). [CrossRef]

**19. **D. F. Walls and G. J. Milburn, 2nd ed.*Quantum Optics* (Springer, Berlin, 2008). [CrossRef]

**20. **D. Slepian, “Some comments on Fourier analysis, uncerlainty and modeling,” SIAM Review **25**, 379–393 (1983). [CrossRef]

**21. **D. L. Donoho and P. B. Stark, “Uncertainty principles and signal recovery,” SIAM J. Appl. Math. **49**, 906–931 (1989). [CrossRef]

**22. **E. J. Candès and C. Fernandez-Granda, “Towards a mathematical theory of super-resolution,” arXiv[1203.5871v1] (2012).