What are the advantages of ghost imaging? Multiplexing for x-ray and electron imaging

TL;DR

This paper analyzes the advantages of multiplexing measurement strategies like ghost imaging, Fourier transform spectroscopy, and Hadamard transform crystallography in x-ray and electron imaging, highlighting when they improve resolution and reduce radiation dose.

Contribution

It provides a quantitative framework for comparing multiplexed measurement schemes and identifies conditions under which multiplexing offers significant benefits.

Findings

Multiplexing can improve resolution and signal-to-noise ratio in specific situations.

Sparse signals enhance the advantages of multiplexing, allowing less radiation use.

The framework enables comparison of error and radiation dose across different measurement schemes.

Abstract

Ghost imaging, Fourier transform spectroscopy, and the newly developed Hadamard transform crystallography are all examples of multiplexing measurement strategies. Multiplexed experiments are performed by measuring multiple points in space, time, or energy simultaneously. This contrasts to the usual method of systematically scanning single points. How do multiplexed measurements work and when they are advantageous? Here we address these questions with a focus on applications involving x-rays or electrons. We present a quantitative framework for analyzing the expected error and radiation dose of different measurement scheme that enables comparison. We conclude that in very specific situations, multiplexing can offer improvements in resolution and signal-to-noise. If the signal has a sparse representation, these advantages become more general and dramatic, and further less radiation can be used to complete a measurement.

Tables

Table 1. Table 1: Quantitative comparison between the sampling schemes discussed, produced using equations ( 4 ) and ( 6 ). Recall m 𝑚 m is the number of measurements, n 𝑛 n is the number of channels (i.e. unknowns), a 𝑎 a is the average dose per channel, σ 𝜎 \sigma is the standard deviation of the additive Gaussian noise per measurement. The factor 1.66 1.66 1.66 in the half Gaussian distribution approximates ( 1 − 2 / π ) − 1 / 2 superscript 1 2 𝜋 1 2 \left(1-2/\pi\right)^{-1/2} . Note that all these values assume A T ​ A superscript 𝐴 𝑇 𝐴 A^{T}A is not singular, which can only occur when m ≥ n 𝑚 𝑛 m\geq n . ‡ Hadamard error estimate derived in [ 14 ] .

	Dose ($m$ meas.)	Error (per channel, $m$ meas.)	Dose-Error Product
Raster	$am$	$\sigma /\left(a\sqrt{mn}\right)$	$\sigma \sqrt{m/n}$
Half Gauss	$\sqrt{2/\pi }\cdot amn$	$1.66\cdot \sigma /\left(a\sqrt{m}\right)$	$1.32\cdot \sigma n\sqrt{m}$
Bernoulli	$\left(1/2\right)\cdot amn$	$2\cdot \sigma /\left(a\sqrt{m}\right)$	$\sigma n\sqrt{m}$
Hadamard‡	$\left(1/2\right)\cdot amn$	$2\cdot \sigma /\left(a\sqrt{m}\right)$	$\sigma n\sqrt{m}$

Equations

$${\mathbf{r}}_m=\sum _i{a}_{mi}\cdot {\mathbf{x}}_i$$

$$R=AX.$$

$$R=AX+ϵ.$$

$$\hat{X}=\underset{X}{\mathrm{argmin}}||R-AX|{|}_2^2,$$

$$\hat{X}=(A^{T}A{)}^{-1}{A}^{T}R,$$

$$\mathcal{E}=\sigma\sqrt{n^{-1}\Big{\langle}\mathrm{Tr}\big{\{}(A^{T}A)^{-1}\big{\}}\Big{\rangle}_{A}}\,.$$

$$\sigma =\sqrt{{⟨{\sigma }_p^2\sum _{i=1}^n{a}_{mi}+{\sigma }_x^2\sum _{i=1}^n{a}_{mi}^2+{\sigma }_m^2⟩}_A}$$

$$\sigma \approx \{\begin{array}{ll}{\sigma }_p\sqrt{n⟨a_{mi}⟩}& \mathrm{Poisson}\\ {\sigma }_x\sqrt{n⟨a_{mi}^2⟩}& \mathrm{Channel}\mathrm{Gaussian}\\ {\sigma }_m& \mathrm{Measurement}\mathrm{Gaussian}\end{array}$$

$$\mathcal{D}={⟨\sum _{ij}|A_{ij}|⟩}_A.$$

$$\mu (A,\Psi )=\sqrt{n}\underset{i,j}{\max }\frac{{\mathbf{a}}_i\cdot {\psi }_j}{||{\mathbf{a}}_i|{|}_2\cdot ||{\psi }_j|{|}_2}$$

$$\underset{X_s}{\min }||R-A\Psi X_s|{|}_2^2+\lambda ||X_s|{|}_1.$$

$$m\ge c{\mu }^2\cdot s\log n$$

$$||\mathbf{r}-A\mathbf{x}|{|}_2^2=$$

$$=$$

$$-2A^T\mathbf{r}+2A^{T}A\mathbf{x}$$

$$\mathbf{r}=A\mathbf{x}+ϵ$$

$$\hat{\mathbf{x}}-\mathbf{x}=A^{+}ϵ$$

$${\mathcal{E}}^2=$$

$$=$$

$$\hat{\bm{\mathbf{x}}}-\langle x\rangle_{i}=\Big{\langle}\bm{\mathbf{a}}_{i}\cdot(r_{i}-\langle r\rangle_{i})\Big{\rangle}_{i}$$

$$\hat{\mathbf{x}}-⟨x{⟩}_i=A^T(\mathbf{r}-⟨r{⟩}_i).$$

$$\hat{\mathbf{x}}=\frac{n{A}^T}{\mathrm{Tr}{A}^{T}A}\mathbf{r}$$

$$c=\underset{c}{\min }||\hat{\mathbf{x}}-\mathbf{x}|{|}_2^2$$

$$ϵ=$$

$$=$$

$$\delta =$$

$$=$$

$$\delta \cdot \delta =(\mathbf{x}\cdot \mathbf{x})\left[(c\lambda -1)\cdot (c\lambda -1)\right]$$

$${\sigma }^2={\sigma }_p^2\sum _i{a}_{ni}+{\sigma }_x^2\sum _i{a}_{ni}^2+{\sigma }_m^2$$

$$r_m=\sum _i{c}_i\mathrm{for}{c}_i\sim \mathrm{Pois}(a_{mi}\cdot x_{mi})$$

$$r_m={\mathbf{a}}_m\cdot {\mathbf{x}}_m+{ϵ}_m$$

Full text

††preprint: APS/123-QED

I Introduction

There has been recent interest in the use of multiplex sensing techniques in applications at the atomic scale [1, 2, 3, 4, 5, 6, 7, 8]. Traditionally, measurements using x-rays or electrons systematically interrogate a system at a single point in space, time, or photon/electron energy at once before moving on to the next measurement. In contrast, multiplex methods measure a combination of space/time/energy points simultaneously. Many repeated multiplexed measurements can convey the same information as a systematic scan if the points combined are chosen wisely. In some cases a multiplex scheme has clear advantages, either fundamental or practical. In other cases, systematic scanning is preferable. The purpose of this paper is to provide a framework for reasoning about which is best in the context of x-ray and electron science.

Through a very general analysis, we are able to conclude that for signals linear in the x-ray or electron intensity:

Multiplex scanning enables the experimenter to swap the resolution – in time, space, or energy – of a final detector for the resolution with which the input beam can be modulated or measured. 2. 2.

We confirm multiplexing can help overcome specific kinds of noise, a well-known result known as Felgett’s advantage. We show, however, this advantage offers no possibility to acquire more information for a given radiation dose as compared to systematic scanning. 3. 3.

Multiplex measurements should yield significant improvements in experiment time, radiation dose used, and signal-to-noise if the sample is known to be sparse in some basis. Multiplexing methods can lead to very general applications of the theory of compressed sensing, which represents an opportunity for imaging at the atomic scale.

Each of these possible advantages will be discussed in detail, with clear applicability conditions. To aid readability, nearly all mathematics has been relegated to the appendices and only results are presented.

To build some intuition for how multiplex sensing techniques work, consider acquiring an x-ray fluorescence image of a sample. We could go about this in the usual way: by rastering a focused x-ray probe over the sample, thereby illuminating one “pixel” at a time, and recording a spectrum at each position. At the end of the experiment, the fluorescence image is simply the set of all acquired spectra, with each spectrum assigned to a spatial location on the image.

Alternatively, we could de-focus the x-ray beam and modulate the transverse intensity distribution with a random mask. This mask might illuminate half of the pixels from our previous raster scan, chosen at random. Then, a measurement acquires the summed spectra from all the illuminated pixels, which by itself contains no spatial information. Imagine, however, we fabricate a large number of such masks and repeat the experiment many times. Every time a specific pixel $i$ is illuminated, it will contribute the $k$ points in its spectrum $\bm{\mathbf{x}}_{i}$ to the spectral sum; for mask $m$, let’s write that measured sum

[TABLE]

where $a_{mi}$ is 1 if the pixel $i$ is illuminated by the mask $m$, and 0 otherwise. For many repeated measurements this can be nicely written in matrix notation

[TABLE]

This form hints that if we build our mask set $A$ in a smart way and take enough measurements, we may be able to easily solve for $X$, the same per-pixel spectra as we acquired in our raster scan.

In the example presented, we specifically designed the masks and therefore knew exactly how to write our matrix $A$, with 1s and 0s for illuminated or shadowed pixels for each mask. The scheme would work just as well, however, if we had random masks we did not choose specifically, so long as we were able to measure and record which pixels were illuminated. This gives us ability to use the randomness inherent in the experiment to our advantage – so long as we can measure that random patterning to high resolution. Figure 1 shows an illustration of all three methods discussed for acquiring a signal: raster scanning, multiplexing with known patterning, and multiplexing using measured patterning. The example envisioned is spatial in nature. The same procedure works just as well, however, in the time or energy (spectral) domains.

When is multiplexing advantageous? In this paper we provide a framework for analyzing this question, then discuss three specific cases where multiplex sensing may be useful. Our focus is on applications involving the use of x-rays and electrons, where measurements are in many ways expensive: they are often flux-limited, cause irreversible damage to the sample, can be difficult to manipulate, and require time at costly facilities. Such practical concerns will color the discussion.

An important note on the use of prior information in imaging experiments is necessary. In general, if one knows something about the signal about to be acquired, for example if it is smooth or always-positive, that fact can be used to improve the measurement process. These priors can be implemented equally well in any sampling scheme, whether rastering or multiplexed. There is one powerful prior, however, for which this equivalence is known to not be true. If the signal is sparse in some representation (i.e. compressible) then multiplex sensing can offer serious advantages over raster scanning [9, 10]. Because of this unique feature, sparsity will be the only prior discussed in this paper. We will first consider sampling schemes that do not make use of sparsity, then show how if the signal is known to be sparse, that information can be readily incorporated into the analysis of different imaging schemes.

Our analysis builds on progress from a variety of fields. One specific implementation of multiplex sensing is classical “ghost imaging” (Fig. 1), which refers to a class of experiments where a randomly patterned wavefront is split into two branches [11]. One branch measures the total transmission through a sample of interest, called the bucket signal. The other branch images the wavefront. A moment of thought reveals this process to be a case of spatial multiplexing, identical to the imaging fluorescence spectroscopy scan discussed previously, only using random beam patterning instead of known masks. In the ghost imaging community, multiplexing with known masks is conversely known as “computational ghost imaging” [12]. Similarly, Fourier transform spectroscopy is also a multiplexing scheme, where many wavelengths are measured simultaneously per measurement [13, 14]. Our presentation generalizes these examples and many others into a single framework. Our treatment is completely classical in nature, ignoring any information gain that might be achieved via performing measurements with entangled particles. Entanglement has been focus of specific studies in ghost imaging – for perspective see [11] – but is not heavily used in current applications.

While the idea to use multiplexing in x-ray and electron experiments has only recently started to recieve significant attention, demonstrations at synchrotrons [1, 3, 2], FELs [7], and electron sources [6] have already been conducted, paving the way for multiplexing to start to impact applied x-ray and electron science.

II Framework

Consider a measurement process where we perform $m$ measurements with $n$ unknown channels, for example the pixels in an image or amplitudes for given frequencies in a spectrometer. The channel values are the unknowns we wish to infer. We measure the weighted sum of all channels; weights may be zero. Crucially, we assume the system response is linear in all aspects, including as a function of beam intensity and the number/character of the channels. Further, we consider the general case where each channel’s output can be described by a $k$-length vector, for example $k=3$ if a pixel gives a simultaneous scalar output for red, green, and blue when measured, and these colors can be recorded separately. Finally, we assume some additive noise in each measurement.

Our goal is choose a weighting scheme that will allow us to infer each individual channel’s $k$ dimensional response. We write this symbolically as follows. Let $R$ be the $m\times k$ matrix of recorded outcomes, $A$ be the $m\times n$ sensing matrix, and $X$ the $n\times k$ signal of interest. Then we consider a measurement process in the presence of uncorrelated, additive noise $\epsilon$ (an $m\times k$ matrix, but where all $k$ pixels are assumed to have equivalent noise) with finite variance $\sigma^{2}$,

[TABLE]

The details of this additive noise will be discussed in a moment.111We have not considered the case where there are errors in the record of the sensing matrix $A$. Such errors are possible when using random $A$ matrices generated naturally as part of the experiment and then measured by some diagnostic, for example as is done in ghost imaging (see Fig. 1). Sensing matrix errors make our estimates of $X$ too small, an effect known as regression dilution. Regression dilution has been studied in detail in the field of statistics and we refer the interested reader to [26]. We will consider a solution to (1) to be the estimate $\hat{X}$ that minimizes the error function

[TABLE]

appropriate for Gaussian errors on the model predictions. For some experiments, other error models may be more appropriate, e.g. Poisson errors in low-flux situations. A brief comment on extending the results here to those more specific cases is given in Appendix E. In the rest of the text we consider Gaussian errors.

The solution to (2) is given by the ordinary least squares solution,

[TABLE]

from which we can compute the total expected error per channel,

[TABLE]

While the ordinary least squares solutions are well known, for completeness in Appendix A we compactly prove (3) and (4). We are motivated to do so by the fact that in the field of ghost imaging it is still common to use sub-optimal estimators, even though least squares has been employed previously [16, 17, 18]. The properties and performance of the traditionally used ghost imaging estimator is discussed in Appendix B and shown numerically in Fig. 2.

Equation (4) naturally separates the error contributions from the noise, given by $\sigma$, and the stability of the inversion of $A$ represented by the trace term $\mathrm{Tr}(A^{T}A)^{-1}$. As we will see, both depend on the choice of $A$, and are therefore determined by the experimental design. The noise factor, $\sigma$, depends also on the noise characteristics of the detection method, which turns out to be important for the comparison between different sampling schemes.

II.1 Noise Model and Transferred Noise

While in general the presented equations are applicable in any case of additive noise, in real experiments the noise can often be broken into three standard types:

Quantum Poisson noise 2. 2.

Per-channel (e.g. per-pixel) Gaussian noise, often caused by the electronics in a detector 3. 3.

Total-detector Gaussian readout noise that is incurred once per measurement.

Considering these three types of noise with variances $\sigma_{p}^{2}$, $\sigma_{x}^{2}$ and $\sigma_{m}^{2}$ respectively, we can approximate the total additive noise $\sigma$ as a function of the magnitudes of the individual contributions and the experimental design,

[TABLE]

where we have assumed each column $\bm{\mathbf{a}}_{m}$ of the matrix $A$ is identically and independently distributed (Appendix C). We have also assumed that we known nothing about the structure of the sample a priori, and therefore assume a uniform signal intensity purpose of estimating Poisson noise (Appendix C.1).

Typically, one source of noise will be dominant in the system. In that case, the other sources will be negligible and we can approximate,

[TABLE]

Note that for the per-channel cases, the noise scales as $\sqrt{n}$, the square root of the number of channels. This should be intuitive, as adding an additional channel brings with it that channel’s noise contribution. It also brings that channel’s contribution to the summed signal, so these factors cancel with the $n^{-1/2}$ factor in equation (4) to make the overall error $\mathcal{E}$ independent of the number of channels.

More interesting is the case where the error is proportional only to the number of measurements. This occurs for example if there is some significant constant readout noise for a detector – then, the number of channels doesn’t matter, only the number of readouts (measurements). In this case, the overall error is reduced by increasing the number of channels and scales as $\mathcal{E}\sim n^{-1/2}$. The advantage comes from each channel contributing additional intensity, and that intensity boosts the signal higher and higher above the constant noise floor. This is well known in the field of spectroscopy, where the $n^{-1/2}$ reduction of the noise as a function of multiple channels is known as the Felgett advantage, one motivation behind Fourier transform and Hadamard spectroscopies. See [14] for a comprehensive study of noise in multiplex spectroscopy.

The question is if this kind of noise is applicable in modern x-ray and electron imaging setups. Current x-ray detectors are able to routinely count individual photons with high quantum efficiency, which effectively eliminates readout noise of this sort [19]. Electron detectors have also progressed into the quantum counting regime recently [20]. The major source of noise in modern setups, therefore, is quantum (Poisson) in nature, and occurs on a per-channel basis. The traditional Felgett advantage does not apply in such situations, and no signal-to-noise advantage is achieved by multiplexing.

II.2 Inversion Stability: Trace Term

Equation (3) assumes that $A^{T}A$ is non-singular and can be inverted, which in general is not true, especially if $A$ is random. In such cases approximate – but often very good – solutions may be found by solving Eq. (2) via a matrix factorization or a minimization algorithm. Reasonable solutions will only be obtained if $A^{T}A$ is high rank, and good design of a sampling scheme should ensure that $A^{T}A$ is non-singular with high probability in the limit of a reasonable number of samples. The schemes we discuss here are all of this variety. In the case of random sensing matrices, however, $A^{T}A$ may only be strictly non-singular in the limit of a large number of random samples, and errors computed using (4) will be lower bounds.

II.3 Radiation Dose

In this model, the total dose on the sample $\mathcal{D}$ for a set of measurements is given by the sum of all the elements of $A$,

[TABLE]

We assume the damage is linear in total dose, and do not consider cases where the damage mechanism changes as a function of fluence, i.e. changes in dose-intensity per unit time or area.

III The Analysis of Some Examples

We consider four sampling schemes for a demonstration comparison: a raster setup, two random multiplexed schemes, and one deterministic multiplexing scheme. For each, we compute the error of a recovered signal as a function of the number of measurements performed and radiation dose incurred. Results are summarized in Table 1.

III.1 Raster Scans

In a raster scan, each pixel is sequentially illuminated with intensity $a$ so that $A\equiv aI$. Usually the number of measurements equals the number of channels ($m=n$). In the presence of noise, it might be productive to repeat measurements. In that case $A$ once again becomes an $m\times n$ matrix where rows of the identity matrix start to repeat.

The solution to the reconstruction problem (Eq. 3) is simply $\hat{X}=R/a$, reflecting the intuition that in raster scanning no reconstruction is necessary. The radiation dose of the experiment is $am$ while the expected error is $(\sigma\sqrt{n})/(a\sqrt{m})$, showing how increasing the beam intensity overcomes background noise at the cost of possibly damaging the sample.

III.2 Multiplex Scans

Any non-raster scan is a multiplex scan. To highlight the differences, we will discuss situations where the sensing matrix $A$ is fairly dense, i.e. about half or more of the elements are non-negligible. Consider the following three concrete multiplex scans:

Half Gaussian. $a_{ij}\sim|\mathcal{N}(0,a^{2})|$. Each element distributed independently and identically (iid). 2. 2.

Bernoulli with $p=1/2$. $a_{ij}=a$ with probability $1/2$ and [math] otherwise, iid. 3. 3.

Hadamard. The Hadamard sequence used in spectroscopy involves constructing a specific $n=2^{z}-1$ length sequence of [math]s and $1$s (with $z$ a positive integer) that forms an $n\times n$ circulant matrix $S_{n}$. This matrix is used as a sensing matrix.222These matrices are closely related to, but are not exactly, the traditionally defined Hadamard matrices. They are formed by generating a Hadamard matrix of order $n+1$, truncating the first row and column – which contain only 1s – and performing the substitions $1\to 0$ and $-1\to+1$. The sequence is deterministic and the circulant property can simplify the experimental implementation. For our discussion, the unframiliar reader need only know that $S_{n}$ is a deterministic matrix where approximately half of the elements are $1$, the rest are [math], and the inverse of $S_{n}$ is given by $2(n+1)^{-1}\cdot(2S_{n}^{T}-J_{n})\neq S^{T}$. By $J_{n}$ we mean an $n\times n$ matrix with all elements $1$ [13, 14].

Note the first two cases are random – continuous and binary, respectively. The last is binary and deterministic.

III.3 Comparison

A quantitative comparison between these four sampling schemes is presented in Table 1, containing analytical results, and Figure 2 which contains numerical results. Some important observations are immediately noticeable:

In cases where radiation damage is a concern, rastering always offers the best error-to-damage tradeoff. The increased damage for adding multiplexing channels is proportional to the number of channels $n$, while the error reduction offered by the Felgett advantage scales as at most $\sqrt{n}$, and at worst negligible, see Eq. (5). 2. 2.

Hadamard sampling can be considered optimal in terms of maximizing Felgett’s advantage [14]. It is interesting, however, that the performance of a random Bernoulli sampling is nearly identical to the carefully designed Hadamard one. 3. 3.

Our results highlight that the performance of all schemes are not the same. The half Gaussian sampling scheme has a higher error for a given dose when compared to the Hadamard or Bernoulli methods, for instance.

IV Experimental Opportunities for Multiplex Sensing

The presented framework gives us the tools for assessing the fundamental capabilities and costs associated with a particular sampling scheme. Often, however, the benefits of multiplexing are practical. The ability to use flexible illumination in either space or time can alleviate engineering challenges, making new experiments possible. Next, we discuss situations where multiplexing can give performance advantages due to either practical or fundamental concerns.

IV.1 Resolution Improvement

In standard, raster-style imaging, the resolution of the measurement is limited by either the probe or detector resolution. For instance, in the x-ray fluorescence imaging application discussed in the introduction, the final resolution of the hyperspectral image obtained will be limited by the x-ray focus size.

Multiplexing enables one to trade the resolution of the probe or detection system for the resolution at which the probe can be patterned or measured. This trade may make it possible to collect data at higher resolution. Alternatively, it may be a simpler or cheaper option than increasing the probe resolution.

As an example, we recently suggested along with our colleagues that by measuring the random temporal fluctuations inherent in the power of an XFEL pulse, pump-probe experiments with time resolution faster than the pulse duration would be possible. A typical XFEL pulse has a duration on the order of 10-100 fs, but its power over time can be readily measured to $\sim 1$ fs resolution using well established diagnostics such as the LCLS’s XTCAV [22]. This should enable x-ray pump-probe measurements to be conducted down to 1 fs resolution. The alternative of decreasing the pump and probe durations – preparing 1 fs x-ray pulses – is possible, but challenging, and requires new hardware to implement [23, 24, 25]. Therefore for some experiments time-domain multiplexing should offer an attractive alternative.

This highlights an interesting aspect of multiplex sensing. In cases where an experimental setup naturally results in random fluctuations in the measurement, it may be possible to use that randomness rather than fight it.

IV.2 Felgett’s Advantage: More Signal per Measurement

Multiplex scans have multiple channels illuminated at once. Thus, for a system that has a fixed flux-per-channel, multiplexed illumination offers a higher overall flux on the sample. In experiments with per-measurement noise, this enables one to achieve a specific signal-to-noise ratio with fewer measurements. Given the large investments in making brighter x-ray and electron sources and the expense associated with operating those sources, the ability to use more of the incident flux is attractive.

It should be emphasized that the experiment time is reduced precisely because the sample is exposed to more incident radiation. This does not necessarily translate into more information acquired for a given radiation dose. Moreover, the multiplex advantage only exists with per-measurement noise (such as $\sigma_{m}$ in Eq. 5), and there is no gain when using photon-counting detectors.

IV.3 Signal-to-Noise per Dose

Can multiplex sensing enable one to obtain more information for a given radiation dose? Many atomic resolution imaging tasks are radiation-damage limited, and the possibility to use multiplex sensing to overcome that limit has been discussed in recent literature [3, 4, 7, 5]. This discussion has revolved around ghost imaging setups, where it has been suggested that damage might be avoided by decreasing the flux on the object branch and increasing the flux into the reference branch. Since the reference branch radiation does not interact with the object, then it follows damage will be reduced.

This argument needs to consider errors in the measurements performed on both branches. The analysis presented so far corresponds to the case of an infinitely bright reference branch with no error on the matrix $A$. In the case where such error is significant, one expects regression dilution [26] and our work only provides a lower bound error estimate. Here, we consider only the limiting case, where the reference branch is made so bright that the error on the measurement of $A$ is negligible.

Now we only need to consider error on the object branch (as written in Eq. 4). Even with an ideal detector, there will at least be Poisson noise due to the fact we image with quanta (photons or electrons).

Even in this best case scenario, our analysis shows there is no fundamental way classical ghost imaging – or any multiplexing scheme – can alleviate radiation damage (in the absence of priors, see section V). Multiplexing with $n$ channels, corresponding to pixels on the reference branch, reduces the error by at most a factor of $\sqrt{n}$ (Eq. 5) which occurs in the presence of per-measurement noise. From a signal-to-noise perspective, this is equivalent to reducing the number of measurements, and therefore the radiation dose, by a factor of $n$ (since error scales with the square root of the number of measurements, $\mathcal{E}\sim m^{-1/2}$). The $n$ additional channels, however, increase the dose by a factor of $n$ (Eq. 6). Thus, the additional dose used to obtain Felgett’s advantage exactly offsets the dose reduction possible by conducting fewer measurements (see Table 1). In the best case scenario, multiplexing does not improve the error-to-dose ratio.

In the case of a perfect detector with only quantum noise, the scaling is even worse – multiplexing does not decrease the error at all, but still incurs the dose penalty. Alone, multiplexing cannot limit radiation damage, only make it worse.

We can understand this intuitively. Consider the limiting case in a ghost imaging setup, where for each ghost image exposure the object branch is so weak that exactly a single photon travels along that branch (zero photons would provide no information, as the bucket would always read zero). That photon is either absorbed by the object or transmitted and detected. An infinitely bright reference branch provides a probability distribution for where the photon interacted with the sample. What is the most informative possible distribution? Clearly, a delta function – telling us exactly where the imaging photon went – would be most useful. This corresponds precisely to a raster scan, where we send the photon in a pre-determined and known direction and measure the result. Ghost imaging imparts at most an equal amount of information per object-branch photon as a raster scan, and in general less.

V Compressed Sensing

Throughout our discussion we have assumed that no prior information is used in the signal reconstruction process. There is, however, one important situation where prior information about the structure of the signal can – in combination with multiplex scanning – dramatically improve the imaging process.

If the signal $X$ to be measured is sparse in some basis, then techniques of compressed sensing can be applied [9, 10]. Here we quickly review compressed sensing theory, in order to show how multiplexed measurements can enable a general application of compressed sensing in x-ray and electron science.

Compressed sensing presumes that our signal $X$ of interest has a sparse representation. Concretely, this means that there exists some orthonormal basis $\Psi$ such that $X_{s}\equiv\Psi^{-1}X$ is $s$-sparse, meaning it has only $s$ non-negligible ($\approx 0$) entries. If $s\ll n$ we may be able to determine the signal $X$ of interest with order $s$ measurements instead of $n$ (Fig. 3). If we can do so, we can obtain the signal with less time, radiation dose, and error.

The basis $\Psi$ can be chosen to maximize the sparsity. For example, if the signal is composed of a waveform with only a few frequencies, choosing $\Psi$ as a DFT matrix will transform the signal into the frequency domain, where it will be sparse. Further, it is absolutely possible for the signal to be sparse in the original domain of interest (for example, an image with a uniform background), in which case $\Psi=I$ is a reasonable choice.

If we knew precisely which elements in $X_{s}$ were sparse, the procedure to use would be straightforward. We could simply truncate $X_{s}$ to the non-negligible components and perform our inference as before with $s$ instead of $n$ unknowns, under the image formation model $R=A\Psi X_{s}$.333We are assuming a complete column of $k$ elements in $X_{s}$ is zero or not for the purposes of counting sparse elements. If one element is non-zero, that entire column is counted as one of the $s$ datapoints we have to infer. If, however, we do not know which components will be sparse, it is essential that each measurement provides some new information about each non-sparse component. This is where multiplexed measurements enter the picture. To ensure that each of our measurements are highly likely to give information about all $s$ components, the measurement scheme $A$ must be constructed in a particular way. It turns out random multiplexing satisfies this requirement for any signal and any sparsity basis, as we will now discuss.

In the theory of compressed sensing, the ability of a sensing scheme to yield information about all the non-zero signal components is quantified by the coherence $\mu$

[TABLE]

which is the correlation between the rows of $A$ and the basis vectors (columns) of $\Psi$. The term coherence should not be confused with the phase relationship between any waveforms, and to distinguish it we will call it CS-coherence.

Given a pair $A,\Psi$, we can replace Eq. 2 with the lasso regression algorithm,

[TABLE]

Here, the $||X_{s}||_{1}$ term enforces signal sparsity and parameter $\lambda$ determined by the noise level $\epsilon$. In general this parameter is unknown, and set using cross-validation or some other estimate based on observed data. If

[TABLE]

for some constant $c$ that depends on $A$, then a central result of compressed sensing theory demonstrates this program recovers the exact signal with probability $1-O(e^{-\gamma m})$ with $\gamma>0$ [10]. If $s$ is small, then often $m\ll n$ measurements give an exact reconstruction. Figure 4 shows this scaling is achieved using the lasso algorithm for a simple sparse reconstruction. This reduction in the number of measurements is the advantage compressed sensing offers for x-ray and electron experiments.

Consider instead trying to use a raster scan ($A=I$) to reveal a naturally sparse signal ($\Psi=I$). Here, $A\Psi=I$ and so $\mu^{2}(A,\Psi)=n$, confirming the expectation that $n$ independent measurements are required to uniquely infer a signal of size $n$.

Since we wish to presume as little as possible about the sample of interest and its sparsity structure, it is natural to ask: are there sensing matrices $A$ that are CS-incoherent with nearly all sparsity representations $\Psi$? It turns out the answer is yes, and randomly multiplexed schemes, such as the half-Gaussian and Bernoulli matrices discussed in this paper, are examples. It can be shown the random nature of such sensing matrices make it almost certain $\mu\ll n$. In other words, signals that are randomly multiplexed are guaranteed to be CS-incoherent. This holds for any sparsity basis $\Psi$ (specifically $\mu^{2}<c\log n$, and often $\mu=$ const.) [9, 10]. This theory shows why multiplexing is so useful in compressed sensing applications – random multiplexing enables us to use any sparsity basis of our choice, thereby maximizing sparsity and the associated advantages of compressed sensing.

As an intuitive example, consider an example where the signal of interest consists of an image with only a few localized features on blank background, but the feature locations are unknown (here, $\Psi=I$). Raster scanning for the features will be slow. In fact for $n$ pixels, $s$ of which are non-blank, we get useful information only a fraction ($s/n$) of the time we perform a raster measurement. Further, we have to scan all the pixels to ensure we don’t miss any feature. Instead, consider conducting measurements where each pixel has a $50\%$ probability of being illuminated in each measurement. On average, we obtain information on fraction $s/2$ of the information-containing pixels per measurement, as compared to the average fraction $s/n$ in the raster scan. Since it is likely $n\gg 2$, this is a more efficient measurement scheme. Further, the probability of missing a given feature goes as $2^{-m}$, which means it is unlikely we miss any feature entirely, even for a modest number of measurements ($<0.1\%$ for $m=10$).

We emphasize that by reducing the number of samples below $n$, compressed sensing is the only multiplexing approach which reduces the total dose needed to achieve a given signal-to-noise ratio.

VI Summary

We have presented a framework for assessing the performance of different sampling schemes. In addition to the common rastering method, where a single point in space, time or energy is measured, we have considered multiplex schemes where sets of points are measured simultaneously. In specific cases, the latter can offer advantages in experiment time, sample use, and measurement resolution:

• •

When per-measurement readout noise is dominant, the Felgett advantage improves the signal-to-noise ratio.

• •

If the source is naturally random, multiplexed measurements using that randomness may be more convenient or yield higher resolution than producing a raster beam.

• •

If the resolution at which the multiplex beam can be measured or controlled exceeds that of a raster beam, multiplexing will improve the resolution of the final reconstruction.

• •

Finally, if the acquired image is known to be sparse in some domain, compressed sensing enables experiments to be conduced with fewer measurements, improving experiment time, radiation dose, and signal-to-noise. Random multiplexing allows one to apply compressed sensing generally, to any signal, with any sparsity basis.

We have proven two facts that may be especially counterintuitive for this budding community. First, the only way multiplexing can be used to mitigate radiation damage or perform more rapid experiments is when it is combined with compressed sensing. If the signal is known to be sparse, the expected number of measurements for the multiplexing schemes scales roughly as $s\log n$, with simulations providing a more accurate final assessment. Second, random multiplex patterns are about as good as the carefully programmed Hadamard sequence – whichever is easier to implement experimentally should be used. Finally, we emphasize an old result that is little known in photon science; there is no Felgett advantage when Poisson noise dominates.

The advantages offered, however, are idiosyncratic to each application. When considering a multiplex sampling scheme, we recommend answering four questions first:

Is the experimental goal limited by radiation dose, number of measurements, or resolution? 2. 2.

Is it possible to program a multiplex scheme or measure some natural or induced randomness in the measurement to implement multiplexing? 3. 3.

What noise is present in your signal of interest? 4. 4.

Is your expected signal known to be sparse in some basis? If so, what is the expected sparsity?

With the answers to these questions in hand, Eq. 4 provides the final error expected for each possible experimental setup.

Large capital investments are currently being made in x-ray and electron facilities and instrumentation. It is our belief that compressed sensing offers an under-utilized opportunity to push the capabilities of these new radiation sources, making current experiments more efficient and opening new experimental possibilities. Since compressed sensing and multiplexing go hand-in-hand, this in turn means we anticipate an increased use of multiplex sensing in x-ray and electron science in the near future. Future work is necessary to develop general ways to apply compressed sensing to the wide variety of imaging currently being performed with electrons and x-rays.

Even in cases where sparsity cannot be leveraged, we believe multiplexing offers significant advantages for very specific x-ray and electron imaging experiments. Our hope is that the work presented here helps clearly identify those advantageous cases and implement them.

Code

Code used to perform simulations and generate figures used in this manuscript is freely available at https://gist.github.com/tjlane/35a84123e3df73cea014ff082a5ad0a8.

Appendix A OLS Proof

For completeness, we quickly sketch proofs of the OLS solution and error analysis. We use vector notation (k=1) for simplicity, but the results are easily extended. To show the optimality of the estimator (3), we wish to minimize (2),

[TABLE]

differentiation with respect to $\bm{\mathbf{x}}$ yields

[TABLE]

which when set equal to [math] shows that (2) is minimized when $\bm{\mathbf{x}}=(A^{T}A)^{-1}A^{T}\bm{\mathbf{r}}$. Note this estimator is defined for non-square $A$; for such matrices $A^{-1}$ nor $(A^{T})^{-1}$ exist. This operator is also known as the Moore-Penrose pseudoinverse and denoted $A^{+}\equiv(A^{T}A)^{-1}A^{T}$, a compact notation we will use immediately. To derive the expected error of this estimator under the error model (1),

[TABLE]

we can compute the residual ($\hat{\bm{\mathbf{x}}}-\bm{\mathbf{x}}$). To do so, multiply by $A^{+}$ and realize $\hat{\bm{\mathbf{x}}}\equiv A^{+}\bm{\mathbf{r}}$ and $A^{+}A=I$ to obtain

[TABLE]

the expected mean squared error is, by definition, the sum of the squared residuals,

[TABLE]

where in the last step we have expanded $A^{+}=(A^{T}A)^{-1}A^{T}$ and used the circular property of the trace ($\mathrm{Tr}ABC=\mathrm{Tr}BCA$) twice.

Appendix B Traditional Ghost Imaging

A quick remark is in order concerning the “traditional” algorithm for ghost imaging, as our discussion reveals how this method works and describes its performance. In the field of ghost imaging, the signal $\bm{\mathbf{x}}$ is reconstructed via correlation with “bucket” detector measurements $r_{i}$ corresponding to many illumination patterns $\bm{\mathbf{a}}_{i}$. Then the image estimate $\hat{\bm{\mathbf{x}}}$ is obtained via correlation

[TABLE]

in matrix notation

[TABLE]

However the forward model is still $\bm{\mathbf{r}}=A\bm{\mathbf{x}}$. Thus, the traditional ghost imaging algorithm implicitly assumes $A^{T}A\approx c_{1}I+c_{2}J$. The interesting observation is that while in general this assumption is not satisfied, we have shown that error-minimizing multiplexing schemes approach this behavior, explaining the success of the traditional ghost imaging algorithm.

That said, the traditional algorithm is inferior in both theory and practice to the ordinary least squares solution given by Eq. (3). The Gauss-Markov theorem proves this so long as the noise is uncorrelated, zero-mean, and has identical variance for all measurement elements.

Note: to obtain properly normalized estimates, like those presented in Fig. 2, one can use

[TABLE]

where we have assumed $\hat{\bm{\mathbf{x}}}$ and $\bm{\mathbf{r}}$ are centered (mean-subtracted). If we let $\hat{\bm{\mathbf{x}}}=cA^{T}\bm{\mathbf{r}}$, this is the form for which

[TABLE]

Proof. Let $\delta\equiv(cA^{T}A-I)\bm{\mathbf{x}}$. Then by definition,

[TABLE]

expand $\bm{\mathbf{x}}=\sum_{i}b_{i}\phi_{i}$ in the basis of the eigenvectors $\phi_{i}$ of $G\equiv A^{T}A$ (a Grammian and therefore orthogonal matrix),

[TABLE]

now

[TABLE]

which is minimized for $c=n/\sum_{i}\lambda_{i}$.

Appendix C Additive Error

The noise model we consider combines noise from three sources: unavoidable Poisson (quantum) noise, per-pixel Gaussian (Johnson type) noise, and finally per-measurement Gaussian (detector readout) noise. Each is shown to be well approximated by a Gaussian distribution. In this appendix we derive the expected variances of each of these types of noise individually; the final noise can be written as the sum of their variances. As a result we obtain that the error on the measurement $r_{m}=\bm{\mathbf{a}}_{m}\cdot\bm{\mathbf{x}}_{m}$ that is a mean-zero normal with variance

[TABLE]

Here, $\sigma_{p}$ gives the magnitude of the Poisson noise relative to the per-pixel $\sigma_{x}$ and per-readout $\sigma_{m}^{2}$ Gaussian variances.

C.1 Poisson Error

Consider summing photons recorded on a set of pixels, with each pixel’s photon count distributed as $c_{i}\sim\mathrm{Pois}(a_{mi}\cdot x_{mi})$,

[TABLE]

then $r_{m}\sim\mathrm{Pois}(\bm{\mathbf{a}}_{m}\cdot\bm{\mathbf{x}}_{m})$, since the sum of iid Poissons is also a Poisson. The mean and variance of this distribution are both $\bm{\mathbf{a}}_{m}\cdot\bm{\mathbf{x}}_{m}$. For large values of this mean, above $\sim 1000$ (photons), a Gaussian approximation is very accurate, and $r_{m}\sim\mathcal{N}(\bm{\mathbf{a}}_{m}\cdot\bm{\mathbf{x}}_{m},\bm{\mathbf{a}}_{m}\cdot\bm{\mathbf{x}}_{m})$. We could equivalently write

[TABLE]

with $\epsilon_{m}\sim\mathcal{N}(0,\bm{\mathbf{a}}_{m}\cdot\bm{\mathbf{x}}_{m})$. We wish, however, to consider the error for a generic measurement system independent of the structure of the sample under study. To do this we assume $\bm{\mathbf{x}}_{m}\sim\sigma_{p}^{2}\bm{\mathbf{1}}$, with $\sigma_{p}^{2}$ a scale factor. Then the Poisson noise can be approximated by $\epsilon_{m}\sim\mathcal{N}(0,\sigma_{p}^{2}\sum_{i}a_{mi})$.

C.2 Per-Pixel Gaussian Error

First, recall that for $\epsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\sigma_{x}^{2})$ that $S\equiv\sum_{i}a_{i}\epsilon_{i}$ is distributed $S\sim\mathcal{N}(0,\sigma_{x}^{2}\sum_{i}a_{i}^{2})$. This means that if we have per-pixel errors $[\bm{\mathbf{\epsilon}}_{m}]_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\sigma_{x}^{2})$ so that

[TABLE]

we can re-write those errors as a per-measurement error

[TABLE]

where $\epsilon_{m}\sim\mathcal{N}(0,\sigma_{x}^{2}\sum_{i}a_{mi}^{2})$.

C.3 Per-Measurement Gaussian Error

The simplest case is a per-measurement error, where

[TABLE]

and $\epsilon_{m}\sim\mathcal{N}(0,\sigma_{m}^{2})$ is a scalar additive error.

Appendix D IID Sensing Matrices

Here we derive the error for sensing matrices with iid components. Let $A=(a_{ij})$ be such a matrix. Per Eq. (4), our task is to evaluate the expected value

[TABLE]

Write the ($n\times n$) Grammian matrix $G=A^{T}A$ for convenience, the elements of which are

[TABLE]

matrix $G$ has diagonal elements $d\equiv m\langle a_{ii}^{2}\rangle$ and off-diagonal elements $c\equiv m\langle a_{ij}\rangle^{2}$. The eigenvalue equation $G\psi=\lambda\psi$ can be written

[TABLE]

for all $n$ degenerate eigenvalues $\lambda$, which correspond to vectors $\psi$ with a single $\psi_{1}$ element and $n-1$ elements $\psi_{2}$. The solution to this system gives $\psi_{2}/\psi_{1}=-1/n$ and more relevantly

[TABLE]

and so

[TABLE]

the terms of this expression are real squared values and therefore positive. This shows that Grammian matrices that approximate the identity matrix, with large values on the diagonal and small off-diagonal matrices, have the lowest error, as expected.

For the half-Gaussian distribution with variance $a^{2}$, $\langle a_{ii}^{2}\rangle=a^{2}$ and $\langle a_{ij}\rangle^{2}=(2/\pi)a^{2}$ so $\left\langle\mathrm{Tr}G^{-1}\right\rangle=n/(ma^{2})\cdot(1-2/\pi)^{-1}$.

For the Bernoulli distribution that gives value $a$ with probability $p$ and [math] otherwise, $\langle a_{ii}^{2}\rangle=pa^{2}$ and $\langle a_{ij}\rangle^{2}=p^{2}a^{2}$ so $\left\langle\mathrm{Tr}G^{-1}\right\rangle=(n/m)\left[a^{2}p(1-p)\right]^{-1}$.

Appendix E Generalized Linear Models

The loss function (2) is appropriate if the errors of the model are Gaussian. In some cases, however, this is not true and superior results can be obtained using a more appropriate model of the noise. One example that will be familiar to those working in x-ray or electron imaging is the low-flux situation, where error is primarily due to Poisson quantum noise, and the error is skewed as a result.

In such situations, it is possible to employ a generalized model to better model the error. Generalized Linear Models (GLMs) are well studied (as in [28]). We only provide a brief overview here for completeness.

Instead of minimizing the least-squares error, as in (2), one can write a more general log-likelihood function

[TABLE]

Here, $g$ is a “link function” that provides a more general model to represent the likelihood of observing data $\bm{\mathbf{r}}$ given input $A\bm{\mathbf{x}}$, and $\gamma$ is a function capturing (all) prior beliefs that weights $\bm{\mathbf{x}}$ according to their believed likelihood before any data are observed.

The desired solution is then the maximum of the log-likelihood

[TABLE]

Equation (2) is a special case where $g$ is the Euclidean distance given by a Gaussian likelihood and no prior beliefs are asserted ($\gamma=\mathrm{const.}$).

Maximum-likelihood theory shows that the asymptotic distribution of $\hat{\bm{\mathbf{x}}}$ is normally distributed around the true value $\bm{\mathbf{x}}^{*}$,

[TABLE]

with covariance given by the inverse of the Fisher Information

[TABLE]

as previously, the expected mean squared error is

[TABLE]

but $\sum_{i}\left(\hat{x_{i}}-x_{i}\right)^{2}$ is the trace of the covariance matrix, so

[TABLE]

as a concrete and relevant example, consider the aforementioned Poisson regression case. Here, the link function is the exponential and the log likelihood becomes,

[TABLE]

where $A_{i}$ is the $i^{\mathrm{th}}$ row of $A$, corresponding to the $i^{\mathrm{th}}$ measurement. Taking the first derivative gives

[TABLE]

which, when set to zero, gives a convex equation that can be solved numerically (via e.g. gradient ascent) to find the maximum likelihood solution. The second derivative is

[TABLE]

One can in principle design an optimal sensing matrix $A$ by minimizing the trace of the inverse of this matrix. If we want to consider such a design before any knowledge of the signal $\bm{\mathbf{X}}$ is know, we might approximate $A_{i}\cdot\bm{\mathbf{x}}\approx\mathrm{const.}$ in which case we obtain the familiar

[TABLE]

Acknowledgements.

TJL would like to thank Anders Nilsson for generously hosting him at Stockholm University during the time this paper was written. This work was supported by the operating contract for the Linac Coherent Light Source at SLAC National Accelerator Laboratory, DoE Contract No. DE-AC02-76SF00515