Spectral Methods for Passive Imaging: Non-asymptotic Performance and Robustness∗

arXiv:1708.04343v3 [cs.IT] 23 Aug 2017

Kiryung Lee† , Felix Krahmer‡ ,

and

Justin Romberg†

Abstract.
We study the problem of passive imaging through convolutive channels. A scene is illuminated with an unknown,
unstructured source, and the measured response is the convolution of this source with multiple channel responses,
each of which is time-limited. Spectral methods based on the commutativity of convolution, first proposed and
analyzed in the 1990s, provide an elegant mathematical framework for attacking this problem. However, these now
classical methods are very sensitive to noise, especially when working from relatively small sample sizes.
In this paper, we show that a linear subspace model on the coefficients of the impulse responses of the channels can
make this problem well-posed. We derive non-asymptotic error bounds for the generic subspace model by analyzing
the spectral gap of the cross-correlation matrix of the channels relative to the perturbation introduced by noise.
Numerical results show that this modified spectral method offers significant improvements over the classical method
and outperforms other competing methods for multichannel blind deconvolution.
Key words. passive imaging, blind deconvolution, perturbation analysis, random matrices
AMS subject classifications. 15B52, 93B30, 94A12

1. Introduction. We give a rigorous analysis of the passive imaging problem. A scene is illuminated by an ambient source that we cannot control or observe. This source signal is always active,
having no discernible “on” or “off” time, and is unstructured. We observe the convolution of this
source with M unknown channel impulse response sequences over a window of time. The goal is to
estimate this ensemble of impulse responses, which in many applications reveals the structure of the
environment being sensed. Problems of this type arise in a wide variety of applications including
as opportunistic channel estimation in underwater acoustics [4, 41–43], seismic interferometry [8],
and passive synthetic aperture imaging [13].
As described fully in Section 2, this is a multichannel blind deconvolution problem, where we
observe the output of a number of linear time-invariant systems all driven by a common source.
We will focus entirely on estimating the responses of these system, and treat the (unknown) source
signal as a supporting actor whose only role is to help us collect information about these channels.
When the channel impulse responses have a finite length K, then a fundamental technique
for performing this estimation, developed in the signal processing literature in the 1990s (see, for
example, [37,47]), is to form a cross-correlation matrix from the channel outputs and then estimate
the channel responses by estimating the null space of this matrix. This method is reviewed in
Section 2.2 below. This classical theory shows that the cross-convolution method is consistent: as
the number of noisy observations we make increases, the channel estimates asymptotically become
aligned with the true underlying impulse responses.
From a finite number of samples, the stability of this process, both in theory and in practice,
depends critically on the spectral gap of this cross-correlation matrix. For even the simplest concrete
instances of this problem, this gap tends to be vanishingly small; a typical example is shown in
Figure 2 (and described in Sections 2.2 and 2.3 below).
∗
The material in this paper was presented in part at the International Conference on Sampling Theory and Applications
(SampTA), Tallinn, Estonia, July 2017 [28].
Funding: This work was supported in part by NSF grants IIS 14-47879, CCF 14-22540 and by DFG grants KR
4512/1-1, KR 4512/2-1.
†
School of ECE, Georgia Institute of Technology, Atlanta, GA (kiryung@ece.gatech.edu, jrom@ece.gatech.edu).
‡
Department of Mathematics, Technical University of Munich, Garching, Germany (felix.krahmer@tum.de).

1

2

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

The main contribution of this paper is to show that if an additional structural constraint on
the channel responses is imposed, then this spectral gap provably widens, stabilizing this channel
estimation procedure. In particular, we constrain the length-K channel responses to be members
of a known D dimensional subspace. Enforcing this constraint requires only a straightforward
modification to the cross-convolution method. Our results show that for a generic D dimensional
subspace (i.e. a subspace chosen at random), the principal angle between the true channel responses
and their estimates decreases as (a) the number of observations in each channel increases, (b) the
signal-to-noise ratio of the observations increases, (c) the number of channels increases, (d) the
model becomes more restrictive, meaning D decreases relative to K.
Our analysis of the subspace constrained cross-convolution method (SCCC) uses as its starting point the classical Davis-Kahan bound on the deviation of the eigenvectors computed from a
perturbed observation of a positive semidefinite matrix. Bounding the size of this perturbation
in terms of the observation noise involves bounding the spectral norms of random matrices with
entries given as coupled high order polynomials of subgaussian random variables. These norms
are written as the suprema of second order chaos processes, for which there are recently developed
concentration results [25,27]. Application of these concentration results involves computing entropy
estimates for various norms. In particular, the entropy estimate for a block norm in Appendix E
is a novel result derived using the polytope approximation and polar duality, which might be of
independent interest.
The Monte Carlo simulation results in Section 4 demonstrate the practical gains that the SCCC
method offers over the classical cross-convolution method. In practice, the estimator produces
accurate results when the number of samples per channel L is a relatively small multiple of K.
In this regime, the estimation error scales (as a function of L) in the same manner as the oracle
solution, where the source is known and the channels are recovered using standard least-squares.
We also demonstrate that SCCC outperforms both the classical method and the recently proposed
approach in [34] for a underwater acoustics simulation with a realistic (non-random) subspace
model.
Related work. As mentioned above, the multichannel blind deconvolution problem was studied
with intense interest in the signal and image processing literature in the 1990s; the methods most
closely related to the work below are described in [18,37,47], and good overviews of general work on
this problem can be found in [35,44]. Many of these algorithms use models on the source signal and
channels, and develop consistency results under different modeling assumptions; see [14, 17, 19, 48]
for representative examples from image processing. To our knowledge, no theoretical results exist
for these algorithms when there are a finite number of noisy samples. More recently, necessary and
sufficient conditions for the generic identifiability of this problem under various geometric priors
have been presented in [32].
A different linearization for the multichannel problem was introduced in [3, 36, 38] and recently
studied thoroughly in [34]. The model presented there is different in that the channels are not
limited in time, a key piece of beneficial structure that our method exploits. These methods also
impose a structural constraint on the source signal, while we view the source signal as unstructured.
Single channel blind deconvolution of signals belonging to low-dimensional subspaces has also
been rigorously studied recently. Identifiability results under various models were studied in [6,
7, 23, 31, 33]. Convex optimization algorithms based on “lifting” were analyzed in [1], followed
by a similar result for a gradient descent algorithm [30]. An alternating minimization algorithm
for blind deconvolution under sparsity models also with subsampling has been analyzed in [29].
While it is possible to extend these methods to the multichannel scenario, unlike the context of
passive imaging, one needs strong geometric priors on both the source and impulse responses.

SPECTRAL METHODS FOR PASSIVE IMAGING

3

This scenario is different from what we consider in this paper. Simultaneous wavelet estimation
and deconvolution of seismic reflection signals [22], auto-calibrated parallel imaging [16], motion
deblurring using multiple images [49] are examples of relevant applications of multichannel blind
deconvolution with geometric priors.
2. Spectral Methods for Blind Deconvolution. In this section, we formulate the FIR multichannel blind deconvolution and describe spectral methods based on the cross convolution.
2.1. Problem formulation. We observe an unknown signal x ∈ CL convolved with multiple
unknown channel impulse responses h1 , . . . , hm ∈ CL with the observations corrupted by additive
noise wm ∈ CL :
(1)

y m = hm ~ x + wm ,

m = 1, . . . , M,

where the convolution ~ is circular1 , i.e.
ym [`] =

L
X

hm [k]x[(k − `) mod L] + wm [`].

k=1

Our goal is to recover the channel responses {hm } from the observations {ym }.
We will assume that the filters have impulse responses of length K; this simply means that the
last L − K entries of each hm are zero. We denotes these non-zero entries using h ∈ CK , with the
relation


hm = S ∗ hm , S = IK 0K,L−K .
The operator S : CL → CK restricts a given vector of length L to its subvector with the first K
elements. The adjoint S ∗ : CK → CL pads L − K zeros to a given vector of length K.
When L ≥ 2K, the nonzero terms in the linear convolution hm ∗ x and the circular convolution
hm ~ x will match, so this model applies to scenarios when we have fixed channels that are being
continuously excited by an unknown input, and we observe a “snapshot” of length L of their
outputs.
2.2. Cross-convolution method. Our method is a modification of the cross-convolution method
introduced 20 years ago in [47]. The core idea is simple: we use the fact that multiple convolutions commute with one another to impose a set of linear constraints that the channel responses
must obey, and then find the (unique up to scale) set of channels responses that obey these linear
constraints. To see how this is done, suppose that the measurements we make are free of noise,
ym = x ~ hm . Then for any pair of channels m, n
(2)

ym ~ hn = x ~ hm ~ hn = yn ~ hm .

Thus the pair of observations yn , ym can be used to construct a set of L constraints on the 2K
variables in the channel coefficients hn , hm .
1

In our model, the source is opportunistic and “always on”, and so the observations in (1) might better be
modeled by a windowed (time-limited) linear convolution. To make a strict correspondence, the source would need
to be periodic, which is an additional structural assumption. Having access to the full circular convolution greatly
simplifies the analysis, and the discrepancy between these two models is marginal when the number of observations
L dominates the length of the impulses responses K. Spectral methods similar to the one presented above that are
explicitly based on the time-limited linear convolution are presented in [37, 47]; in practice, small gains might be
realized by using these closely related methods.

4

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

To make this more precise, let Tv be the L × K matrix whose action Tv h circularly convolves
v = [v1 , . . . , vL ]> ∈ CL with h after zero-padding:
Tv = C v S ∗ ,
where Cv ∈ CL×L is a circulant matrix defined by


v1
vL vL−1 · · · v2
 v2
v1
vL · · · v3 


Cv :=  .
.
.
.
.

.
.
vL vL−1 vL−2 · · · v1
Then we can write (2) as
Tym hn − Tyn hm = 0.
We can represent all M (M − 1)/2 such constraints in one linear system. With




0L,K . . . 0L,K
Tyi+1
−Tyi
Y (1)
 .

..
..
..
 Y (2) 
 ..

.
.
.


(i)
,
(3) Y = 
 , where Y = 
..


−Tyi
T yM


L,K . . . 0L,K
.
 0

{z
}
|
{z
}
|
Y (M −1)
(M − i) × (i − 1) blocks
(M − i) × (M − i) block diagonal
> >
we know that h = [h>
1 , . . . , hM ] will be in the null space of Y . Indeed, under the mild condition
that the z-transforms of the hm do not share common zeros, the null space of Y is one dimensional,
containing only the scalar multiples of h [47].
In any practical scenario, noise (and possibly other perturbations), will keep (2) from holding
exactly, and Y will in general not have a null space. The channel estimates, then, are formed by
finding the vector that is as close to a null vector as possible; after forming Y from the observations,
we solve

(4)

minimize kY vk22
v∈CM K

subject to kvk2 = 1.

The solution to the above is of course given by the eigenvector of Y ∗ Y corresponding to the smallest
eigenvalue.
The matrix Y can be unwieldy for large M , its dimensions are M (M − 1)/2 × KM . However,
we can form the smaller KM × KM matrix Y ∗ Y in a computationally efficient way using fast
convolutions. Y ∗ Y can be thought of as an M × M array of K × K matrices; a quick calculation
shows that K ×K block Bn,m , corresponding to rows (n−1)K +1 to nK and columns (m−1)K +1
to mK in Y ∗ Y , is given by
(P
∗
m0 6=m Tym0 Tym m = n,
Bn,m =
−Ty∗m Tyn
m 6= n.
Thus Y ∗ Y can be computed with M (M + 1)/2 convolutions of length L. Computing the solution
to (4) can be done with an eigenvalue decomposition in O(M 3 K 3 ) time. For large values of M K,
the solution can be computed with the power method, with each application of Y ∗ Y computed
using fast convolutions.
Under certain statistical assumptions on the noise, this estimate is consistent: as L → ∞, the
smallest eigenvector of Y ∗ Y goes to (a scalar multiple of) h. However, to date there is no rigorous

SPECTRAL METHODS FOR PASSIVE IMAGING

5

analysis of the stability of this procedure. There are no non-asymptotic accuracy bounds that tell
us what kind of performance we should expect for a certain number of channels M , filter lengths
K, and observation times L.
The effect of noise on the accuracy of the estimate given by (4) can be understood using the
spectral properties of the “noise-free” cross correlation matrix. We write the noisy measurements
as
ym = sm + wm , where sm = hm ~ x.
The cross correlation matrix Y is simply the sum of the cross correlation matrix Ys for the signals
sm (i.e. create Ys as in (3) using Tsi in place of the Tyi ) and the cross correlation matrix Yw for
the noise signals wm . The estimate of the channels is formed by solving
(5)

minimize v ∗ (Ys∗ Ys + E) v
v∈CM K

subject to kvk2 = 1,

where
E = Ys∗ Yn + Yn∗ Ys + Yn∗ Yn .
From the discussion above, we know that in the noise-free case (E = 0), we will recover the true
channel responses. In expectation, the E matrix becomes a scalar multiple of the identity, and
the eigenvectors (and relative order of the eigenvalues) does not change. From a finite number of
samples, how closely the solution to (5) matches the noise-free solution depends on the size of E
relative to the spectral gap of Ys∗ Ys , which is the size of its second smallest (or smallest non-zero)
eigenvalue. This is codified in the classical Davis-Kahan sin θ-theorem [10].
Theorem 2.1 (sin θ theorem [15, Corollary 7.2.6]). Let A, E ∈ Cn×n satisfy that A and A + E
b denote the eigenvector of A (resp. A + E) corresponding
are positive semidefinite. Let q (resp. q)
to the smallest eigenvalue. Suppose that λn−1 (A) > λn (A). If
(6)

kEk ≤

λn−1 (A) − λn (A)
,
5

then
(7)

b ≤
sin ∠(q, q)

4kEqk2
.
λn−1 (A) − λn (A)

Remark 2.1. The error bound in the Davis-Kahan theorem is known to be sharp for general
perturbations. Recent results in [39, 46] have provided refined bounds for unstructured random
perturbations, but unfortunately do not apply to our perturbation matrix E.
Since the eigenvectors are unit norm, having a bound on angle between them is almost the
same as having an error bound (up to a global phase), i.e.
√
e ≤ min kq − eiθ qk
e 2 ≤ 2 sin ∠(q, q).
e
(8)
sin ∠(q, q)
θ∈[0,2π)

As discussed above, when the channels are identifiable, λKM (Ys∗ Ys ) = 0, and so we will have
guarantees for the robustness of (5) when λKM −1 (Ys∗ Ys ) is large compared to kEk. Unfortunately,
this smallest non-zero eigenvalue λKM −1 (Ys∗ Ys ) is typically very small in magnitude. Figure 1a
shows a typical example; here we create Ys from M = 4 channels of length K = 256; the channel impulse responses themselves were generated at random. In this example, λ1 (Ys∗ Ys ) = 1 but
λKM −1 (Ys∗ Ys ) = 4.7 × 10−5 , and so we only have robustness guarantees for the mildest perturbations. The practical performance of the estimator is poor in even mild amount of additive noise,
as the experiments in Figure 2 suggest.

6

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

2.3. Subspace-constrained cross-convolution (SCCC) method. In this paper, we show that
the introduction of a linear model for the channel responses can tangibly increase the size of this
spectral gap. Using a linear subspace to model the channel responses has had some empirical success
in the literature. For example, in [43] a data-driven linear model is constructed for underwater
acoustic channels for the purpose of ocean tomography.
Along with having an impulse response of limited length, we will also assume that the hm
lie in known subspaces of dimension D < K. This means that each hm can be expressed as
hm = Φm um , where the columns of Φm form a basis for the model subspace, and the um are the
expansion coefficients in this basis — recovering the um is now the same as recovering the channel
responses hm . The concatenated channels are written

h = Φu,

(9)


Φ1

Φ=


..


,

.
ΦM




u1


u =  ...  .
uM

With this model in place, the channel coefficients u will be in the null space of Ys Φ.
The estimation procedure has to be modified to account for a slight bias introduced by the
∗ ] = σ 2 I, we have
linear model. With random uncorrelated noise, E[wm ] = 0, E[wm wm
w
E[Yn ] = 0,

2
and E[Yn∗ Yn ] = σw
(M − 1)LI,

and so
2
E[Φ∗ Y ∗ Y Φ] = Φ∗ Ys∗ Ys Φ + σw
(M − 1)L · Φ∗ Φ.

To make the perturbation from the noise-free cross-correlation matrix zero mean, we will solve
(10)

2
minimize z ∗ Φ∗ (Y ∗ Y − σw
(M − 1)LIM K )Φz,
z∈RM D

subject to kzk2 = 1.

2 (M −
Again, the solution is the eigenvector corresponding to the smallest eigenvalue of Φ∗ (Y ∗ Y −σw
1)LIM K )Φ.

Remark 2.2. If Φ∗ Φ = IM D , then since adding a scalar multiple of the identity does not perturb
2 (M − 1)LI
eigenvectors, we may ignore σw
M K in (10). Otherwise, subtracting the noise covariance
2
∗
σw (M − 1)LIM K from Y Y further suppresses the error in the estimated impulse responses. In
2 needs to be estimated and the error in this estimate will propagate
practice, the noise variance σw
2 is known a
to the estimate of the impulse responses. For simplicity of analysis, we assume that σw
priori.
Figure 1b shows the effect of the subspace constraint on the spectral gap. Here, a generic
subspace was chosen by generating Φ at random. The entries of Φ ∈ C256×8 were generated as
independent copies of a standard complex Gaussian random variable. The size of smallest non-zero
eigenvalue is now significantly more distinct (λM D−1 (Φ∗ Ys∗ Ys Φ)/λ1 (Φ∗ Ys∗ Ys Φ) = 0.4). As the
numerical results in Section 4 show, adding subspace constraints of this nature does indeed lead to
significant robustness of the method in the presence of noise.
Our main results, detailed in Section 3, quantify this spectral gap for generic subspaces Φ.
3. Main Results.

SPECTRAL METHODS FOR PASSIVE IMAGING

7

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0
200

400

600

800

1000

(a) Ys∗ Ys

10

20

30

(b) Φ∗ Ys∗ Ys Φ

Figure 1. (a) The eigenvalue spectrum of an example cross correlation matrix created from M = 4 channels of
length K = 256 with random coefficients. Note that there is almost no gap between the second smallest eigenvalue
and zero. (b) The eigenvalue spectrum after introducing a linear model of dimension D = 30. Note that the spectral
gap is now pronounced.

3.1. Non-asymptotic Analysis. Our main results give non-asymptotic performance guarantees
for the subspace-constrained cross-convolution method. We make the following two assumptions
throughout:
(A1) Generic subspaces. It is clear that some choices of linear channel models will be better
than others. We will investigate the “generic” case, where the bases themselves are generated at random. In particular, we will assume that Φ1 , . . . , ΦM are independent copies of a
K-by-D complex Gaussian matrix whose entries are independent and identically distributed
(iid) as CN (0, 1). Our theorems below hold with high probability with respect to this draw
of the Φm ; we might interpret this as saying that the results hold for “most” subspace
models. In Section 4 below, we empirically confirm this performance also for choices of Φm
with additional structure as it would appear in applications, even though they do not arise
from the model analyzed in this paper.
(A2) Random noise. The perturbations to the measurements w1 , . . . , wM ∈ CL are iid sub∗ ] = σ 2 I , and are independent of the bases
gaussian vectors with E[wm ] = 0 and E[wm wm
w L
{Φm }.
We present two main theorems below. In the first one, we assume that the input itself is a
white random process. In the second one, we study deterministic inputs with a mild structural
assumption on the common source signal that essentially amounts to x being spread out in the
frequency domain; the resulting error bounds are slightly weaker than for the random model.
The theorems provide sufficient conditions on the number of samples L we need to observe at
the output of each channel in order to guarantee a certain level of accuracy in the estimate ĥ found
b and then taking ĥ = Φb
by solving (10) for u
u. The number of samples we need will depend on the
length of the filter responses K, their intrinsic dimensions D, the number of channels M , and the
signal-to-noise-ratio (SNR) defined as
(11)

P
2
Eφ [ M
m=1 khm ~ xk2 ]
η :=
.
P
2
Ew [ M
m=1 kwm k2 ]

Under (A1) and (A2), it follows from the commutativity of convolution and Lemma B.1 that η
simplifies as
(12)

η=

Kkxk22 kuk22
.
2
M Lσw

8

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

In addition, the bounds will depend on the channel impulse responses all being roughly the same
size. We measure the disparity in impulse response energies using the flatness parameter
√
M kum k2
(13)
.
µ := max
1≤m≤M
kuk2
When a small
√ number of the M impulse responses are significantly greater than the others, we will
have µ ≈ M . In this case, we expect to have longer observation times, as we are only getting a
small number of diverse looks at the signal. Our results are most interesting when µ is a constant
on the order of 1. Qualitatively, this means that each channel is roughly as important as the others.
We now present the first of our main results. Theorem 3.1 below assumes a random common
source signal x. The bound on the number of observations L sufficient to guarantee a certain
accuracy in the channel estimates is a complicated expression involving the number of channels M ,
their maximum impulse response lengths K, their intrinsic dimension D, the SNR η, the channel
flatness µ, and the level of accuracy . But in reasonable scenarios where the noise is not extreme
(η is a constant), and K and M are not too different, we have
√
b h) ≤  when L & KD/,
sin ∠(h,
with the inequality on the right holding to within log factors. As a point of reference, we are
estimating M D channel coefficients from M L samples at the outputs; we have the same number
observations as unknowns when L & D. As D ≤ K, our estimate scales in a mildly unsatisfying
way, though as the recovery procedure is highly nonlinear, it is unclear what form an optimal
scaling would take.
Theorem 3.1 (Random Source). We observe noisy channel outputs {ym } as in (1), with SNR η
as in (11), and form an estimate h of the channel responses by solving (10). Suppose assumptions
(A1) and (A2) above hold, let x be a sequence of zero-mean iid subgaussian random variables with
variance σx2 , η ≥ 1, µ = O(1), and L ≥ 3K.2 Then for any β ∈ N, there exist absolute constants
C > 0, α ∈ N and constants C1 (β), C2 (β) such that if there are a sufficient number of channels,
(14)

M ≥ C1 (β) logα (M KL),

that are sufficiently long,
(15)

K ≥ C1 (β)D logα (M KL),

and we have observed the a sufficient number of samples at the output of each channel,

C1 (β) logα (M KL)  K
(16)
L≥
+
D
,
η
M2
then with probability exceeding 1 − CK −β , we can bound the approximation error as
√
 1  √K √ 
D 
b h) ≤ C2 (β) logα (M KL) √
(17)
sin ∠(h,
+ D + √
.
ηL M
η ML
Remark 3.1. The SNR requirement η ≥ 1 was introduced to simplify the expressions in Theorem 3.1. The conditions in the low SNR regime η < 1 can be easily extracted from the proof of the
theorem and Proposition 3.3 below.
Without the subspace prior, L > K is necessary to claim that Y ∗ Y has nullity 1 in the noiseless case. We used
L ≥ 3K in the proof in order to use the identity that the circular convolutions of three vectors of length K modulo
L indeed coincide with their linear convolution.
2

SPECTRAL METHODS FOR PASSIVE IMAGING

9

Theorem 3.1 is interpreted as follows: When the dimension D of the subspaces in (9) is small
(up to a fraction of the ambient dimension K), the number of channels M is large (depending
weakly on the other dimension
parameters only through log), and the length of observation is large
√
enough (L ≥ 3K + C KD for an absolute constant C), we can apply the Davis-Kahan theorem
b The error bound in (17) converges to 0 if either L or η
which provides an error bound for h.
grows toward infinity. Moreover, the error bound is nonasymptotic since it explicitly shows how
the error depends on L and η when they are finite. In a heuristic argument that counts the number
of unknown parameters and the number of given equations, a necessary condition for the unique
recovery of x and {hm }M
m=1 from noise-free measurements is given as L ≥ M D/(M − 1). It is
unclear whether this is also a valid necessary condition under the finite impulse response structure.
On the other hand, with some diversity in x (not necessarily under an explicit stochastic model), it
has been shown L ≥ 3K suffices in this scenario. Still, it is not clear whether this remains a valid
necessary condition with the extra subspace model in (9). However, empirically, spectral methods
including the classical and what we propose in this paper break down when L < K. Along the
above discussion, it is still open to answer whether the requirement on L in (16) is near optimal or
not.
To prove Theorem 3.1, we establish an intermediate result for the case where the input signal
x is deterministic. In this case, our bounds depend on the spectral norm ρx of the (appropriately
restricted) autocorrelation matrix of x,
e ∗ Cx S
e∗ k,
ρx := kSC
x
where


0K−1,L−K+1
IK−1
e
S=
.
I2K−1
02K−1,L−2K+1

(18)

Then the deterministic version of our recovery result is:
Theorem 3.2 (Deterministic Source). Suppose that the same assumptions hold as in Theorem 3.1, only with x as a fixed sequence of numbers obeying
ρx ≤ C3 kxk22 .

(19)
If (15) and (14) hold, and
(20)

L≥


C1 (β) logα (M KL)  K 2
+
KD
,
η
M2

then with probability exceeding 1 − CK −β , we can bound the approximation error as
(21)


α
(M KL)  K √
b h) ≤ C2 (β) log
√
sin ∠(h,
+ KD .
M
ηL

The condition (19) can be interpreted as a kind of incoherence condition on the input signal x.
Since
ρx ≤ kCx k2 = Lkb
xk2∞ ,
b ∈ CL is the normalized discrete Fourier transform of x, it is sufficient that x̂ is approxiwhere x
mately flat for (19) to hold. This is a milder assumption than imposing an explicit stochastic model
on x as in Theorem 3.1. For the price of this relaxed condition, the requirement on L in (20) that
activates Theorem 3.2 is more stringent compared to the analogous condition (16) in Theorem 3.1.

10

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Theorems 3.1 and 3.2 distinguish from a recent result [34] in the following sense. Ling and
Strohmer [34] analyzed the error bound for the least squares solution to a different linearized
formulation in [3]. In their analysis, the unknown filters were assumed to follow stochastic subspace
models, which span vectors fully supported on the entire observation period. Obviously, these
models do not explain the FIR structures arising in applications. Unlike their analysis, we explicitly
considered the case where the unknown filters are supported on a short interval. On the other hand,
the number of observations L enabling the error bound in [34] scales near optimally whereas L grows
faster in Theorems 3.1 and 3.2. Again, models considered in these analyses are different and it is still
open to verify whether a near optimal scaling can be achieved with limited randomness satisfying
the FIR structure.
3.2. Proof of Main Results. The main results in Theorems 3.1 and 3.2 are obtained by the
following proposition. Proposition 3.3 identifies a sufficient condition for (6), which enables Theorem 2.1 and provides an error estimate of computing the most dominant eigenvector from a noisy
matrix. The sufficient condition is stated in terms of scaling of key parameters for the sake of
interpretation.
Proposition
3.3. Suppose the assumptions in (A1) and (A2) hold, ρx satisfies (19), L ≥ 3K,
√
and µ ≤ M /2.3 Let ρx,w denote the cross-correlation among the input x and the noise terms
w1 , . . . , wM defined by
e ∗ Cwm S
e∗ k,
ρx,w := max kSC
x
1≤m≤M

e ∈ R(3K−2)×L is as in (18). For any β ∈ N, there exist absolute constants C > 0, α ∈ N
where S
and constants C1 (β), C2 (β) that only depend on β, for which the following holds. If
(22)

K ≥ C1 (β)µ4 D logα (M KL),

(23)

M ≥ C1 (β)µ4 logα (M KL),

and
(24)

  D2
 D
ρ2x,w
L
K 
2
µ
≥
+
+
D
+ 2
2 kxk2
C1 (β) logα (M KL)
K
M2
η
ηKσw
2

then
b h)
sin ∠(h,
(25)


≤ C2 (β) log (M KL) √
α

√ 
√
  D
√ 
ρx,w
K
D 
µ √ +
+ D + √
M
ηKLσw kxk2
K
η ML

holds with probability 1 − CK −β .
b = Φb
b of u. Then h
Proof of Proposition 3.3. Recall that we first compute an estimate u
u serves
as an estimate of h = Φu. Since the estimation error is measured in the principal angle, which
is invariant under scalar multiplication, without loss of generality, we may assume that kuk2 = 1.
b and u
b are related by
Indeed, the errors in the estimates h
(26)
3

b =
sin ∠(h, h)

Pb⊥ h
h

khk2

2

≤

b 2
b k2
kh − hk
σmax (Φ)ku − u
σmax (Φ) √
b ),
≤
≤
· 2 sin ∠(u, u
khk2
σmin (Φ)kuk2
σmin (Φ)

By definition, the parameter µ always satisfies 1 ≤ µ ≤

√
√
M . In this perspective, µ ≤ M /2 is a mild condition.

SPECTRAL METHODS FOR PASSIVE IMAGING

11

where the last step follows from (8).
By the assumption in (A1) and (22), the standard results on singular values of subgaussian
matrices (e.g., see [9, Theorem II.13]) imply that the condition number of Φm is upper bounded
by 3 for m = 1, . . . , M with high probability as we choose C1 (β) in (22) large enough. We proceed
the proof under this event. Then the condition number of Φ is also upper bounded by 3.
b in the principal angle.
Therefore, it suffices to focus on bounding the estimation error in u
e = Φ∗ (Y ∗ Y − σ 2 (M − 1)LIM K )Φ.
b is computed as the least dominant eigenvector of A
Note that u
w
2 (M −1)LI
Furthermore, the target vector u is the unique null vector of A = E[Φ∗ (Y ∗ Y −σw
M K )Φ].
To see this, we decompose Y as Y = Ys + Yn , where the noise-free portion Ys (resp. the noise
portion Yn ) is obtained as we replace ym = hm ~ x + wm in Y by its first summand hm ~ x (resp.
by its second summand wm ) for all m = 1, . . . , M . Consequently, we have
2
Ew [Yn∗ Yn ] = σw
(M − 1)LIM K

and Ew [Ys∗ Yn ] = 0

as well as
e = Ew [Φ∗ (Y ∗ Y − σ 2 (M − 1)LIM K )Φ]
Ew [A]
w
= Φ∗ Ys∗ Ys Φ + Φ∗ Ew [Ys∗ Yn ]Φ + Φ∗ Ew [Yn∗ Ys ]Φ
2
+ Ew [Φ∗ (Yn∗ Yn − σw
(M − 1)LIM K )Φ]

= Φ∗ Ys∗ Ys Φ.
As shown in Section 2.2, by the construction of Ys , the vector h = Φu with the true filter
coefficients is in the null space of Ys . Therefore, u is almost surely a null vector of the noise-free
e = E[Φ∗ Y ∗ Ys Φ]. The uniqueness
matrix Φ∗ Ys∗ Ys Φ and hence also of its expectation A = E[A]
s
follows from the first part of the following lemma, which is proved in Section 5.1.
Lemma 3.4. Under the hypothesis of Proposition 3.3, the following are true: i) The nullity of
E[Φ∗ Ys∗ Ys Φ] is 1; ii) Nonzero eigenvalues of E[Φ∗ Ys∗ Ys Φ] are no less than K 2 kxk22 kuk22 /2 =: δ.
This lemma also establishes a lower bound for the gap between the two smallest eigenvalues of A.
This spectral gap allows to distinguish the corresponding eigenspaces of A. Provided condition (6),
e does not deviate too much from its expectation A in the spectral norm (this will be the
that is, A
e
main task of the remainder of this proof), this property also carries over to the eigenspaces of A
e and A are close to each
and it follows from Theorem 2.1 that the least dominant eigenspace of A
b is a good estimate of u.
other. Thus, up to a global phase, u
It remains to show that condition (6) is satisfied with high probability. To this end, we derive
e − A and show that the
a tail estimate of the spectral norm of the random perturbation E = A
perturbation relative to the spectral gap δ satisfies
r 
hr 1
kEk
D
α
≤ C(β) log (M KL)
+
µ2
δ
M
K
√ 
√ i
(27)
  D
√ 
ρx,w
K
D
+√
µ √ +
+ D + √
M
ηKLσw kxk2
K
η L
with probability 1 − CK −β , where C(β) is a constant depending only on β. By choosing C1 (β) in
(22), (23), and (24) large enough, we can make the right hand side of (27) less than 1/5. Thus (6)
is satisfied.
The derivation of (27) is rather involved for the following reasons: The entries of the perturbation matrix E are given as fourth order polynomials of subgaussian random variables. In addition,

12

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

the convolution structure in the construction of Y creates dependence relations between the matrix entries. To analyze the perturbation, we decompose E into three components of different
polynomial order as follows.
E = Φ∗ Ys∗ Ys Φ − E[Φ∗ Ys∗ Ys Φ]
|
{z
}
∗

+Φ
|
+

Es
∗
Ys Yn Φ +Φ∗ Yn∗ Ys Φ

{z

Ec
∗
Φ (Yn∗ Yn

|

}
2
− σw
(M − 1)LIM K )Φ .
{z
}
En

The following lemmas, the proofs of which will be presented in Section 5, provide tail estimates
of the components; the tail estimate in (27) is then obtained by combining these results via the
triangle inequality.
Lemma 3.5. Suppose that (A1) holds. For any β ∈ N, there exist a numerical constant α ∈ N
and a constant C(β) that depends only on β such that
r 
r 1
kΦ∗ Ys∗ Ys Φ − E[Φ∗ Ys∗ Ys Φ]k
D
α
(28)
≤ C(β) log (M KL)
+
µ2
2
2
2
M
K
K kxk2 kuk2
holds with probability 1 − CK −β .
Lemma 3.6. Suppose that (A1) holds. For any β ∈ N, there exists a constant C(β) that depends
only on β such that, conditional on the noise vector w,
√ 
√ 
C(β)ρx,w   D
kΦ∗ Ys∗ Yn Φk
K
√
√
(29)
+
D
≤
µ
+
M
K 2 kxk22 kuk22
ηKLσw kxk2
K
holds with probability 1 − CK −β .
Lemma 3.7. Suppose that (A1) holds. For any β ∈ N, there is a constant C(β) that depends
only on β such that
r
2 (M − 1)LI
kΦ∗ (Yn∗ Yn − σw
D
C(β) logα (M KL)
M K )Φk
(30)
≤
·
2
2
2
η
L
K kxk2 kuk2
with probability 1 − CK −β .
Finally, under the event where (6) is satisfied, Theorem 2.1 implies that
(31)

b) ≤
sin(u, u

4kEuk2
.
δkuk2

To estimate the right hand side, we again decompose E as in (28), so the triangle inequality
yields
(32)

kEuk2 ≤ kEs uk2 + kEc uk2 + kEc∗ uk2 + kEn uk2 .

To bound the first term, recall that u is in the null space of Φ∗ Ys∗ Ys Φ, so we obtain that
(33)

Es u = 0.

SPECTRAL METHODS FOR PASSIVE IMAGING

13

For the second and third summand, Lemma 3.6 yields that with probability 1 − CK −β
√ 
√ 
C(β)ρx,w K 3/2 kxk2 kuk32   D
K
∗
√
µ √ +
+ D .
(34) max(kEc uk2 , kEc uk2 ) ≤ kEc kkuk2 ≤
M
ηLσw
K
A bound for the last summand is provided by the following lemma, which is proved in Section 5.5.
Lemma 3.8. Suppose that (A1) holds. For any β ∈ N, there is a constant C(β) that depends
only on β such that
r
2 (M − 1)LI
kΦ∗ (Yn∗ Yn − σw
D
C(β) logα (M KL)
M K )Φuk2
(35)
≤
·
2
4
2
η
ML
K kxk2 kuk2
with probability 1 − CK −β .
Inserting the bounds for the four summands into (32) yields the error bound in (25), which
completes the proof.
In the remainder of this section, we show how Theorems 3.1 and 3.2 can be deduced from
Proposition 3.3.
Proof of Theorem 3.2. Since most assumptions of the Theorem agree with the ones of Proposition 3.3 it only remains to bound ρx,w . This is achieved by the following lemma, which is proved
in Appendix C.1.
Lemma 3.9. Suppose (A2) holds and let x be a fixed sequence of numbers obeying (19). For any
β ∈ N, there exists an absolute constant C such that
√ p
ρx,w ≤ CKσw ρx 1 + log M + β log K
holds with probability 1 − K −β .
The theorem follows from a direct application of Proposition 3.3.
Proof of Theorem 3.1. We again need to show a bound for ρx,w , but in addition we need to
estimate ρx , as it is not part of the assumptions. The following lemma, which is proved in Appendix C.2, provides these two bounds.
Lemma 3.10. Suppose (A2) holds and let x be a sequence of zero-mean iid subgaussian random
variables with variance σx2 . Then
√
L + Cβ KL log5 (KL)
ρx
√
≤
kxk22
L − 2Lβ log K
and

√
Cβ KL log5 (M KL)
ρx,w
≤ p
√
σw kxk2
L − 2Lβ log K

hold with probability 1 − 3K −β .
Again, the theorem follows from a direct application of Proposition 3.3.
4. Numerical Results. We compare the classical cross-convolution (CC) method and our modification with additional subspace prior, which is the subspace-constrained cross-convolution (SCCC)
method in a set of Monte Carlo simulations.
Our first tests concern the random signal model of Theorem 3.1 with additional subspace constraints. As expected, our method SCCC, which exploits this information, significantly outperforms

14

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

10 0

10

10 0

CC
SCCC
NBLS

-1

10

10 0

CC
SCCC
NBLS

-1

10

-1

CC
SCCC
NBLS

10 -2

10 -2

10 -3

10 -2

10 -3
0.05

0.1

0.15

(a)

0.2

0.25

10 -3
10

20

(b)

30

40

2

4

6

8

10

(c)

b h) for cross-convolution (CC) and subspace-constrained cross-convolution
Figure 2. Estimation error sin ∠(h,
(SCCC). (95th percentile for 1,000 trials). I.i.d. Gaussian basis. Default parameters: K = 256, M = 4, D = 8,
L = 20K, SNR = 20 dB. (a) For different dimensions. (b) For different observation lengths. (c) For different
numbers of channels.

the original CC, which does not, see Figure 2. Specifically, the estimation error monotonically decreases (resp. increases) with L and M (resp. D). This is consistent with the prediction in
Theorem 3.1.
Figure 3 compares the empirical phase transition of SCCC and the least squares approach in
the non-blind
p case (where x is known). As in the limit when M goes to infinity, (16) simplifies
to L/K & D/K, we provide the plot in terms of the quantities D/K and L/K. In the case
of noisy measurements, our performance measure is the 95th percentile of the estimation error,
i.e., we consider the worst case, but ignore those 5% of the instances, which performed worst. In
Figure 3, we compare this error to the non-blind case. Our color coding uses a logarithmic scale
with blue denoting the smallest and red the largest estimation error within the regime. We observe
that, unlike the non-blind case, SCCC totally fails when D/K is larger than certain threshold
determined by M . This phenomenon clearly justifies the need to introduce a strong subspace prior
to stabilize the eigenvector estimation.
Finally we study data obtained from a parametric channel impulse response model and apply
SCCC under subspace model obtained empirically by principal component analysis [43]. More
precisely, the unknown filters are generated by sampling a known continuous function with random
shifts (not necessarily on a given grid) followed by scaling with random amplitudes. We compare
our method to the classical cross-convolution (CC) method as well as to the least squares (LS)
approach to a different linearized formulation [3], which also incorporate the same prior models
on the impulse responses. As shown in Figure 4, SCCC outperforms CC and LS in this scenario.
Although the assumptions of Theorem 3.1 are not satisfied, similarly to the previous experiment,
the estimation error for SCCC monotonically decreases with L and M . The other two methods did
not perform satisfactorily even under a very high SNR of 80 dB. We have already explained why the
classical method fails in terms of the spectral gap. For the least squares methods, which recovers

SPECTRAL METHODS FOR PASSIVE IMAGING

40

0

15

40

40

30

30

20

20

10

10

-0.5

30
-1

20
-1.5

10
-2

0

-2.5

0

0.25

0.5

0.75

0

1

0
0

0.25

(a)

0.5

0.75

1

0

(b)
40

40

30

30

30

20

20

20

10

10

10

0

0
0.25

0.5

(d)

0.75

1

0.5

0.75

1

0.75

1

(c)

40

0

0.25

0
0

0.25

0.5

(e)

0.75

1

0

0.25

0.5

(f)

Figure 3. Estimation error for SCCC (upper row) and nonblind deconvolution with known x (bottom row). (95th
percentile of the log of the estimation error for 1,000 trials). I.i.d. Gaussian basis. x-axis: D/K. y-axis: L/K.
K = 256, SNR = 20 dB. (a), (d) M = 2. (b), (e) M = 4. (c), (f ) M = 6.

both the input and filters simultaneously, was not successful because the system of convolution
with multiple channels is highly ill-conditioned. Even when the unknown filters are known, the
corresponding system has condition number typically larger than 5,000. This happened since the
known continuous function is close to a strict band-pass filter and the unknown signal has a white
spectrum. Figure 4c demonstrates that even under moderate SNRs, SCCC provides stable recovery
whereas the other methods totally failed in this regime.
5. Proof of Key Lemmas. In this section, we prove some important lemmas required in the
proofs of our main theorems. For this it will be of particular importance to compute tail estimates
of the spectral norms of several structured random matrices with entries given as functions of second
order polynomials of Gaussian variables. In certain cases, such bounds can be established via matrix
concentration inequalities (e.g., [21, 24, 45]). These matrix concentration inequalities, however,
require upper estimate of spectral norms of summands and high order moments, the computation
of which turns out rather complicated for those matrices arising in the proofs. Moreover, there are
cases where these inequalities do not apply as the target matrix is not decomposed as a sum of
independent variables. For these reasons, we rather proceed by writing the spectral norms of these
random matrices in a variational form as suprema of corresponding chaos processes. Tail estimates
for such suprema of second order chaos processes, as given in the following theorem, have been

16

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

10 -1

10 -1

10 0

10 -1
10

-2

10

-2

10

10 -3

-2

10 -3
10 -3
CC
SCCC
LS

10 -4

10 -4
10

20

30

40

20

(a)

40

60

10 -4
40

50

(b)

60

70

80

(c)

b h) for cross-convolution (CC) and subspace-constrained cross-convolution
Figure 4. Estimation error sin ∠(h,
(SCCC) (95th percentile for 1,000 trials). PCA basis. Default parameters: K = 64, M = 64, D = 6, L = 20K, SNR
= 80 dB. (a) For different observation lengths. (b) For different numbers of channels. (c) For different SNR.

derived using chaining arguments [25].
Theorem 5.1 (Tail estimates for suprema of chaos processes [25, Theorem 3.1]). Let ξ ∈ Cn be an
L-subgaussian vector with E[ξξ ∗ ] = In , ∆ ⊂ Cm×n , and 0 < ζ < 1. There exists a constant C(L)
that only depends on L such that
sup |kM ξk22 − E[kM ξk22 ]| ≤ C(L)(K1 + K2

p
log(2ζ −1 ) + K3 log(2ζ −1 ))

M ∈∆

holds with probability 1 − ζ, where K1 , K2 , and K3 are given by



K1 := γ2 (∆, k · k) γ2 (∆, k · k) + dF (∆) + dF (∆)dS (∆),
K2 := dS (∆)(γ2 (∆, k · k) + dF (∆)),
K3 := d2S (∆).
Using the polarization identity, this result on the suprema of second order chaos processes has
been extended from a subgaussian quadratic form to a subgaussian bilinear form [27].
Theorem 5.2 (A corollary of [27, Theorem 2.3]). Let ξ ∈ Cn be an L-subgaussian vector with
E[ξξ ∗ ] = In , ∆2 , ∆1 ⊂ Cm×n , 0 < ζ < 1, and a > 0. There exists a constants C(L) that only
depends on L such that
sup
M1 ∈∆1 ,M2 ∈∆2

e1 + K
e2
hM1 ξ, M2 ξi − E[hM1 ξ, M2 ξi] ≤ C(L)(K

p

e 3 log(8ζ −1 )),
log(8ζ −1 ) + K

SPECTRAL METHODS FOR PASSIVE IMAGING

17

e1, K
e 2 , and K
e 3 are given by
holds with probability 1 − ζ, where K



e 1 := aγ2 (∆1 , k · k) + a−1 γ2 (∆2 , k · k) aγ2 (∆1 , k · k) + a−1 γ2 (∆2 , k · k) + adF (∆1 ) + a−1 dF (∆2 )
K



+ adF (∆1 ) + a−1 dF (∆2 ) adS (∆1 ) + a−1 dS (∆2 ) ,



−1
−1
−1
e
K2 := adS (∆1 ) + a dS (∆2 ) aγ2 (∆1 , k · k) + a γ2 (∆2 , k · k) + adF (∆1 ) + a dF (∆2 ) ,

2
e 3 := adS (∆1 ) + a−1 dS (∆2 ) .
K
A special case of Theorem 5.2 where a = 1 was shown in [27, Theorem 2.3]. Note that the
bilinear form satisfies
hM1 ξ, M2 ξi = haM1 ξ, a−1 M2 ξi, ∀a > 0.
Moreover, the γ2 functional and the radii with respect to the Frobenius and spectral norms are all
1-homogeneous functions. Therefore, Theorem 5.2 is a direct consequence of [27, Theorem 2.3].
Since a > 0 in Theorem 5.2 is arbitrary, one can minimize the tail estimate over a > 0.
5.1. Proof of Lemma 3.4. Note that Φ∗ Ys∗ Ys Φ ∈ CM D×M D is an M -by-M block matrix,
where the size of each block is D-by-D. Then it follows from the special structure of Y (with each
row consisting only of some Tyk in position j and the corresponding −Tyj in position k) that the
(m, m)th diagonal block of Φ∗ Ys∗ Ys Φ is given by
(36)

M
X

∗
Tx~hm0 Φm
Φ∗m Tx~h
m0

M
X

=

m0 =1
m0 6=m

m0 =1
m0 6=m

e ∗ C∗
Φ
m Φ
e

m0 u m0

Cx∗ Cx CΦ
e

m0 u m0

e m,
Φ

e m = S ∗ Φm for m = 1, . . . , M . Similarly, the (m, m0 )th off-diagonal block of Φ∗ Y ∗ Ys Φ for
where Φ
s
m 6= m0 is given by
(37)

∗
e ∗ C∗
Tx~hm Φm0 = −Φ
− Φ∗m Tx~h
m Φ
e
m0

m0 um0

e
Cx∗ Cx CΦ
e m um Φm0 .

By Lemma B.3, the expectation of the (m, m)th diagonal block is given by
M
h X
e ∗ C∗
E
Φ
m Φ
e
m0 =1
m0 6=m

m0 u m0

Cx∗ Cx CΦ
e

m0 u m

M
i
X
em =
Φ
kum0 k22 K 2 kxk22 ID
0
m0 =1
m0 6=m

and a similar calculation yields that the expectation of the (m, m0 )th diagonal block for m0 6= m is
given by
e ∗ C∗
−E[Φ
m Φ
e

m0 u m0

2
2
∗
e
Cx∗ Cx CΦ
e m um Φm0 ] = −K kxk2 um um0 .

Collecting the above expectations, we obtain that E[Φ∗ Ys∗ Ys Φ] satisfies
E[Φ∗ Ys∗ Ys Φ]
= kuk22 Pu⊥ − Υ,
K 2 kxk22
where
Υ=

M
X
m=1

em e∗m ⊗ kum k22 Pu⊥
,
m

18

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

where ⊗ denotes the Kronecker product and e1 , . . . , eM denote the standard basis vectors in RM .
Since Υu = 0 and hence u is in the kernel of E[Φ∗ Ys∗ Ys Φ], the gap between its two smallest
eigenvalues is given by its smallest non-zero eigenvalue, which, using that kΥk ≤ max1≤m≤M kum k22 ,
is estimated from below by

max1≤m≤M kum k22  K 2 kxk22 kuk22
K 2 kxk22 kuk22 1 −
≥
.
2
kuk22

(38)

Here the last inequality follows from our assumption that
√
√
M kum k2
M
≤
µ = max
.
1≤m≤M
kuk2
2
> >
5.2. Proof of Lemma 3.5. We use the following notation in the proof. For u = [u>
1 , . . . , uM ] ,
D
where um ∈ C for k = 1, . . . , M , we define a block (p, q) norm by

kukp,q =

(39)

 M
X
1/q



kum kqp

1 ≤ q < ∞,

m=1



 max kum kp
1≤m≤M

q = ∞.

By (36) and (37), Φ∗ Ys∗ Ys Φ−E[Φ∗ Ys∗ Ys Φ] is rewritten as the sum of its diagonal block portion
(a) and the off-diagonal block portion (b), where
(40)
M
M

 X
X
∗
∗
∗
e
e
e∗ ∗
e ∗ C∗
(a) =
em em ⊗
Φ
e 0 u 0 Φm ] ,
e 0 u 0 Φm − E[Φm CΦ
m Φ
e u Cx Cx CΦ
e u Cx Cx CΦ
m=1

(b) = −

M X
M
X
m=1 m0 =1
m0 6=m

m0 =1
m0 6=m

m0

e ∗ C∗
em e∗m0 ⊗ (Φ
m Φ
e

m0

m

m0 u m0

m0

m

e
e∗ ∗
Cx∗ Cx CΦ
e m um Φm0 − E[Φm CΦ
e

m0

m0 um0

m

e
Cx∗ Cx CΦ
e m um Φm0 ]).

Therefore, by the triangle inequality, we have
kΦ∗ Ys∗ Ys Φ − E[Φ∗ Ys∗ Ys Φ]k ≤ k(a)k + k(b)k.
It remains to compute tail estimates for the spectral norms of (a) and (b).
Diagonal block portion: The (m, m)th block of (a) is written as Φ∗m Υs,m Φm , where
Υs,m =

M
X
m0 =1
m0 6=m

∗
SCΦ
e

m0 um0

Cx∗ Cx CΦ
e

m0 um0

m

S∗.

Then due to the block diagonal structure, we have
k(a)k = max kΦ∗m Υs,m Φm − E[Φ∗m Υs,m Φm ]k.
1≤m≤M

Since Φm and Υs,m are independent, E[Φ∗m Υs,m Φm ] is rewritten as
E[Φ∗m Υs,m Φm ] = EΦm [Φ∗m E{Φm0 }m0 6=m [Υs,m ]Φm ] = E{Φm0 }m0 6=m [tr(Υs,m )ID ].

SPECTRAL METHODS FOR PASSIVE IMAGING

19

Therefore, we have
(41)

Φ∗m Υs,m Φm − E[Φ∗m Υs,m Φm ] = Φ∗m Υs,m Φm − tr(Υs,m )ID + (tr(Υs,m ) − E[tr(Υs,m )])ID .
|
{z
} |
{z
}
(c)

(d)

We will compute tail estimates of the spectral norms of (c) and (d) separately and then combine
them using the triangle inequality. First, we compute a tail bound of k(c)k using the following
lemmas, which are proved in Appendices D.1 and D.2.
Lemma 5.3. Let Ψ = [ψ1 , . . . , ψD ] ∈ CK×D satisfy that vec(Ψ) follows CN (0KD,1 , IKD ), where
> ]> . Then
vec(Ψ) = [ψ1> , . . . , ψD
√
kΨ∗ AΨ − E[Ψ∗ AΨ]k ≤ CkAk KD log(8ζ −1 )
holds with probability 1 − ζ.
Lemma 5.4. Suppose that (A1) holds. For any β ∈ N, there exist a numerical constant α ∈ N
and a constant C(β) that depends only on β such that
max kΥs,m − E[Υs,m ]k ≤ C(β)ρx Kkuk2,∞ kuk2 logα (M KL)

1≤m≤M

holds with probability 1 − K −β .
Lemma 5.3 is a direct consequence of the theory of suprema of second order chaos processes
[25, 27]. By Lemma 5.3, conditioned on Υs,m ,
√
k(c)k ≤ C1 kΥs,m k KD(log M + β log K)
holds with probability 1 − M −1 K −β . Then by Lemmas 5.4 and B.1 with the triangle inequality, it
follows that
√
k(c)k ≤ C(β)(ρx + kxk22 )K 3/2 Dkuk22 log2 (M K).
holds with probability 1 − M −1 K −β .
Next we consider (d). Note that
∗
tr(SCΦ
e

m0 u m0

Cx∗ Cx CΦ
e

m0 um0

S ∗ ) = KkCx S ∗ Φm0 um0 k22 .

Therefore, the spectral norm of (d) is rewritten as
M
X

k(d)k =

KkCx S ∗ Φm0 um0 k22 − E[KkCx S ∗ Φm0 um0 k22 ]

m0 =1
m0 6=m

=K

M
X

∗
2
>
∗
2
k(u>
m0 ⊗ Cx S )φm0 k2 − E[k(um0 ⊗ Cx S )φm0 k2 ]

m0 =1
m0 6=m

=K

M
X

∗
(em0 e∗m0 ⊗ u>
m 0 ⊗ Cx S ) φ

m0 =1
m0 6=m

|

{z

Am

}

2
2

h
−E

M
X

∗
(em0 e∗m0 ⊗ u>
m0 ⊗ Cx S )φ

m0 =1
m0 6=m

2i
2

,

20

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

> >
where φm = vec(Φm ) for m = 1, . . . , M and φ = [φ>
1 , . . . , φM ] . In fact, we are computing a
2
tail bound of the Gaussian quadratic form kAm φk2 , which can be done by the Hanson-Wright
inequality, Lemma A.2. Observing that the block diagonal matrix U with blocks um satisfies
kU ∗ U kF ≤ kuk22,4 , we obtain

k(d)k = K|kAm φk22 − E[kAm φk22 ]|
≤ C3 K(kA∗m Am kF ∨ kAm k2 )(2 log M + 2β log K)
√
≤ C3 K( Kkuk22,4 ∨ kuk22,∞ )kSCx∗ Cx S ∗ k(2 log M + 2β log K)
≤ C4 K 3/2 kuk22,4 ρx (2 log M + 2β log K)
holds with probability 1 − M −1 K −β . Note that the tail bound of (c) dominates that of (d).
Collecting the above estimates, it follows that
√
k(a)k ≤ C(β)ρx K 3/2 Dkuk22 logα (M KL)
holds with probability 1 − K −β .
If we normalize with the spectral gap given in (38), then the relative perturbation due to (a) is
upper bounded by
r
k(a)k
D
α
(42)
≤ C(β) log (M KL)
2
2
2
K
K kxk2 kuk2
with probability 1 − K −β .
Off-diagonal block portion: Unlike the diagonal block portion (a), the off-diagonal block portion
(b) does not have a block diagonal structure and computing its tail bound is more involved.
To restrict the convolution of two short vectors of length K to its support, we introduce S̆ ∈
(2K−1)×L
R
defined by


0K−1,L−K+1
IK−1
S̆ =
.
IK
0K,L−K
Then we have
S S̆ ∗ S̆ = S

(43)
and
∗
CΦ
e

(44)

m um

e m = S̆ ∗ S̆C ∗
Φ
e
Φ

m um

e m,
Φ

∀m = 1, . . . , M.

Due to the commutativity of product of two circulant matrices, (43), and (44), we can rewrite
(b) as
(b) = −

M X
M
X



em e∗m0 ⊗ (Cx S̆ ∗ Zm )∗ Cx S̆ ∗ Zm0 − E[(Cx S̆ ∗ Zm )∗ Cx S̆ ∗ Zm0 ] ,

m=1 m0 =1
m0 6=m

where
∗
Zm = S̆CΦ
e

m um

e m,
Φ

m = 1, . . . , M.

Note that the summation in (b) runs over all distinct pairs (m, m0 ) with m 6= m0 . Our main trick
here is to add and subtract the terms corresponding to pairs (m, m) for m = 1, . . . , M . This ends

SPECTRAL METHODS FOR PASSIVE IMAGING

21

up with a diagonal sum and a full summation over all pairs (m, m0 ). The resulting full summation
term now provides a nice factorization, which leads to an analysis using the techniques for the
second-order chaos processes.
Indeed, since the Φm ’s are independent, we have
E[Cx S̆ ∗ Zm ]∗ E[Cx S̆ ∗ Zm0 ] − E[(Cx S̆ ∗ Zm )∗ Cx S̆ ∗ Zm0 ] = 0D,D ,

∀m 6= m0 .

Therefore, (b) is decomposed as (b) = (e) − (f), where
(e) =

(f) =

M
X

em e∗m ⊗ ((Cx S̆ ∗ Zm )∗ Cx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ]∗ E[Cx S̆ ∗ Zm ]),

m=1
M
X

em e∗m0 ⊗ ((Cx S̆ ∗ Zm )∗ Cx S̆ ∗ Zm0 − E[Cx S̆ ∗ Zm ]∗ E[Cx S̆ ∗ Zm0 ]).

m,m0 =1

Note that (Cx S̆ ∗ Zm )∗ Cx S̆ ∗ Zm0 − E[Cx S̆ ∗ Zm ]∗ E[Cx S̆ ∗ Zm0 ] is decomposed as
(Cx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ])∗ (Cx S̆ ∗ Zm0 − E[Cx S̆ ∗ Zm0 ])
∗
+ E[Cx S̆ ∗ Zm
](Cx S̆ ∗ Zm0 − E[Cx S̆ ∗ Zm0 ])

(45)

+ (Cx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ])∗ E[Cx S̆ ∗ Zm0 ].
Therefore, the spectral norm of (e), which corresponds to the extra diagonal term, is upper-bounded
by
k(e)k ≤



+2

max kCx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ]k

2

1≤m≤M





max kCx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ]k
max kE[Cx S̆ ∗ Zm ]k .

1≤m≤M

1≤m≤M

By Lemma B.2, we have
kE[Cx S̆ ∗ Zm ]k = KkCx e1 u∗m k ≤ Kkxk2 kum k2 ≤ Kkxk2 kuk2,∞ .
We again use bounds for suprema of second order chaos processes [25, 27] to get a tail bound
for kZm − E[Zm ]k, as given in the following lemma, which is proved in Appendix D.3.
Lemma 5.5. Suppose that (A1) holds. For any β ∈ N, there exist a numerical constant α ∈ N
and a constant C(β) that depends only on β such that
sup kZm − E[Zm ]k ≤ C(β)kuk2,∞ K logα (M KL)
1≤m≤M

holds with probability 1 − K −β .
By Lemma 5.5, it follows that the relative perturbation due to (e) is upper bounded by
(46)
with probability 1 − K −β .

C(β) logα (M KL)µ2
k(e)k
≤
M
K 2 kxk22 kuk22

22

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Similarly, (f), which corresponds to the full 2D summation, is rewritten as
M
X

(f) =

e∗m

∗

⊗ Cx S̆ Zm

M
X

e∗m ⊗ E[Cx S̆ ∗ Zm ]



=

∗

⊗ Cx S̆ Z

e∗m0 ⊗ E[Cx S̆ ∗ Zm0 ]

e∗m

+

∗

m0 =1
M
∗  X

∗
∗
∗
∗
∗
0
0
em ⊗ E[Cx S̆ Zm ]
em0 ⊗ (Cx S̆ Zm − E[Cx S̆ Zm ])
m0 =1

m=1

+

M
X



M
∗  X

⊗ (Cx S̆ Zm − E[Cx S̆ Zm ])
e∗m0 ⊗ (Cx S̆ ∗ Zm0 − E[Cx S̆ ∗ Zm0 ])
∗

m=1
M
X



m0

m0 =1

m=1
M
X

e∗m0

m0 =1
M
∗  X

m=1

−

M
∗  X

M
∗  X

e∗m ⊗ (Cx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ])
e∗m0 ⊗ E[Cx S̆ ∗ Zm0 ] .
m0 =1

m=1

Therefore, by the triangle inequality, we have
k(f)k ≤

+2

M
X

e∗m ⊗ (Cx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ])

m=1
M
X

2

M
X

e∗m ⊗ (Cx S̆ ∗ Zm − E[Cx S̆ ∗ Zm ])

m=1

e∗m ⊗ E[Cx S̆ ∗ Zm ] .

m=1

> ]> . Then, by Lemma B.2, we have
Let v1 , . . . , vM ∈ CD and v = [v1> , . . . , vM
M
X

e∗m

∗

⊗ E[Cx S̆ Zm ] = sup
v∈B2M D

m=1

M
X
m=1

= sup K
v∈B2M D

E[Cx S̆ ∗ Zm vm ]

M
X
m=1

≤ sup Kkxk2
v∈B2M D

Cx e1 u∗m vm
M
X

2

2

= sup Kkxk2
v∈B2M D

M
X

u∗m vm

m=1

kum k2 kvm k2 ≤ Kkxk2 kuk2 .

m=1

On the other hand,
M
X
m=1

e∗m

∗

∗

⊗ (Cx S̆ Zm − E[Cx S̆ Zm ]) = Cx S̆

∗

M
X


e∗m ⊗ (Zm − E[Zm ])

m=1

≤ kS̆Cx∗ Cx S̆ ∗ k1/2

M
X

| m=1

e∗m ⊗ (Zm − E[Zm ]) .
{z

(‡)

}

√
As implied by (19), the first factor is bounded by C3 kxk2 . Hence it remains to show an upper
bound on the last term (‡). This is established in Lemma 5.6 using the results on suprema of
second-order chaos processes [25, 27] together with an entropy bound by polytope approximation
of a unit ball [20], polar duality, and entropy duality [2]; see Appendix D.4 for the proof.

SPECTRAL METHODS FOR PASSIVE IMAGING

23

Lemma 5.6. Suppose that (A1) holds. For any β ∈ N, there exist a numerical constant α ∈ N
and a constant C(β) that depends only on β such that
M
X

e∗m ⊗ (Zm − E[Zm ]) ≤ C(β)kuk2,∞ (K +

√

M KD) logα (M KL)

m=1

holds with probability 1 − K −β .
Collecting the above estimates and noting that, up to log factors, both factors are bounded by
µKkxk2 kuk2 , we obtain that the relative perturbation due to (f) is upper bounded by
r 
r 1
k(f)k
D
α
≤ C(β) log (M KL)
(47)
+
µ2 ,
2
2
2
M
K
K kxk2 kuk2
potentially for an increased value of α, with probability 1 − CK −β .
Finally, (28) follows by combining (42), (46), and (47).
5.3. Proof of Lemma 3.6. The proof of Lemma 3.6 is similar to (and easier than) that of
Lemma 3.5. We will reuse some of tail estimates obtained in Section 5.2. On the other hand,
we also need to derive tail estimates of the suprema of certain Gaussian processes, which did not
arise in Section 5.2. We use a moment-version of Dudley’s inequality [12] to compute these tail
estimates.
Similarly to the previous section, we use the following decomposition into the diagonal block
portion and the off-diagonal block portion:
Φ∗ Ys∗ Yn Φ = (g) + (h),
where
(g) =

M
X

em e∗m

⊗

M
 X
m0 =1
m0 6=m

m=1

(48)
(h) = −

M X
M
X
m=1 m0 =1
m0 6=m

e ∗ C∗
Φ
m Φ
e

m0 u m0

e ∗ C∗
em e∗m0 ⊗ Φ
m Φ
e


em ,
Cx∗ Cwm0 Φ

m0 u m0

e m0 .
Cx∗ Cwm Φ

We derive upper bounds on the spectral norms of (g) and (h) respectively in the following.
Diagonal block portion: Note that (g) is a block diagonal matrix and its expectation is 0. Define
(49)

Υc,m =

M
X
m0 =1
m0 6=m

∗
SCΦ
e

m0 u m0

Cx∗ Cwm0 S ∗ .

Then it follows from the block diagonal structure that
k(g)k ≤ max kΦ∗m Υc,m Φm k.
1≤m≤M

Since Φm and Υc,m are independent, Φ∗m Υc,m Φm is rewritten as
(50)

Φ∗m Υc,m Φm = Φ∗m Υc,m Φm − tr(Υc,m )ID + tr(Υc,m )ID .
|
{z
} |
{z
}
(i)

(j)

24

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

First, we compute a tail estimate of k(i)k. Similarly to the previous section, we use Lemma 5.3
conditioned on Υc,m . Then we apply the tail estimate of kΥc,m k given in following lemma, whose
proof is in Appendix D.5.
Lemma 5.7. Suppose that (A1) holds. For any β ∈ N, there exist a numerical constant α ∈ N
and a constant C(β) that depends only on β such that, conditional on the noise vector w,
√
max kΥc,m k ≤ C(β)ρx,w Kkuk2 logα (M KL)
1≤m≤M

holds with probability 1 − K −β .
By Lemmas 5.3 and 5.7, we obtain that
√
k(i)k ≤ C(β)ρx,w K Dkuk2 logα (M KL)
holds with probability 1 − M −1 K −β .
On the other hand, a direct calculation shows that (j) is expressed as
tr(Υc,m ) =

M
X

∗
∗
(u>
m0 ⊗ Kwm0 Cx S )vec(Φm0 ),

m0 =1
m0 6=m

where tr(Υc,m ) denotes the complex conjugate of tr(Υc,m ). The above expression implies that
tr(Υc,m ) is a linear function of the independent Gaussian matrix entries, and hence is a zero-mean
Gaussian random variable. The variance of tr(Υc,m ) is given by
M
X

∗ 2
2
2 2
∗
K 2 kum0 k22 kwm
0 Cx S k2 ≤ K kuk2 ρx,w .

m0 =1
m0 6=m

Therefore, by a tail estimate of a Gaussian variable,
p
k(j)k ≤ Cρx,w Kkuk2 1 + log M + β log K
with probability 1 − M −1 K −β . Hence k(i)k dominates k(j)k.
By collecting the above estimates, we obtain that the relative perturbation due to (g) is upper
bounded by
r
√
C(β)ρx,w logα (M KL) D
ρx,w
k(g)k
C(β) logα (M KL)
D
√
√
(51)
≤
=
2
2
2
2
K
K kxk2 kuk2
kxk2 kuk2
ηL
Kσw kxk2 M
with probability 1 − CK −β , where in the last step we used (12).
Off-diagonal portion: Similarly to the analogous part of the proof of Lemma 3.5, we add and
subtract the diagonal sum and obtain
(h) = (k) + (l),
where
(k) =
(52)

M
X

em e∗m ⊗ Φ∗m SCx∗ Cwm S̆ ∗ Zm ,

m=1

(l) = −

M
X
m,m0 =1

em e∗m0 ⊗ Φ∗m SCx∗ Cwm S̆ ∗ Zm0 .

SPECTRAL METHODS FOR PASSIVE IMAGING

25

By the triangle inequality,
k(k)k ≤

M
X

em e∗m ⊗ Φ∗m SCx∗ Cwm S̆ ∗ (Zm − E[Zm ])

|m=1
M
X

+

{z

}

(m)

em e∗m ⊗ Φ∗m SCx∗ Cwm S̆ ∗ E[Zm ] .

|m=1

{z

}

(n)

We use the result by Davidson and Szarek [9, Theorem II.13] to get a tail estimate of kΦm k.
Specifically, it follows from (A1) that
p
√
√
√
√
(53)
max kΦm k ≤ k[ 2Re(Φm ), 2Im(Φm )]k ≤ K + 2D + 2 log M + 2β log K
1≤m≤M

holds with probability 1 − K −β .
By Lemma 5.5 and (53),
k(m)k = max kΦ∗m SCx∗ Cwm S̆ ∗ (Zm − E[Zm ])k
1≤m≤M




≤
max kSCx∗ Cwm S̆k
max kΦm k
max kZm − E[Zm ]k
1≤m≤M

1≤m≤M

1≤m≤M

≤ ρx,w C(β)K 3/2 kuk2,∞ log6 (M K)
holds with probability 1−K −β , where k·k2,∞ is defined in (39). On the other hand, by Lemmas B.2
and A.3,
k(n)k = K max kΦ∗m SCx∗ wm u∗m k
1≤m≤M

≤ Kkuk2,∞ max kΦ∗m SCx∗ wm k2
1≤m≤M

∗
= Kkuk2,∞ max k(wm
Cx S ∗ ⊗ ID )vec(Φ>
m )k2
1≤m≤M
p
∗
≤ Kkuk2,∞ max kwm
Cx S ∗ ⊗ ID kF 1 + log M + β log K
1≤m≤M
p
√
∗
= Kkuk2,∞ max
Dkwm
Cx S ∗ k2 1 + log M + β log K
1≤m≤M
√ p
≤ 2Kkuk2,∞ ρx,w D 1 + log M + β log K
∗ C S∗k ≤
holds with probability 1 − K −β , where the last inequality follows from the fact that kwm
x
2
∗
∗
kSwm Cx S k. Note that k(m)k dominates k(n)k.
The spectral norm of (l) is upper bounded through a factorization by

k(l)k ≤

M
X

e∗m

⊗

∗
S̆Cw
Cx S ∗ Φm
m

{z

}

(o)
M
X

+

e∗m0 ⊗ (Zm0 − E[Zm0 ])

m0 =1

m=1

|

M
X

{z

|

(p)

em e∗m0 ⊗ Φ∗m SCx∗ Cwm S̆ ∗ E[Zm0 ] .

m,m0 =1

|

{z

(q)

}

}

26

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

The spectral norm of (o) is written as the supremum of a Gaussian process and is bounded by
the following lemma, which is proved in Appendix D.6.
Lemma 5.8. Suppose that (A1) holds. For any β ∈ N, there exists a constant C(β) that depends
only on β such that, conditional on the noise vector w,
M
X

∗
e∗m ⊗ S̆Cw
Cx S ∗ Φm ≤ C
m

p
√
√
1 + βρx,w ( M D + K) log K

m=1

holds with probability 1 − K −β .
By Lemma 5.6,
k(p)k ≤ C(β)kuk2,∞ (K +

√
M KD) logα (M KL)

with probability 1 − K −β , where k · k2,∞ is defined in (39). Note that k(o)kk(p)k dominates k(m)k.
Therefore, we may ignore k(m)k.
By Lemma B.2, the spectral norm of (q) is upper bounded by
2

k(q)k =

M
X

em e∗m0 ⊗ KΦ∗m SCx∗ wm u∗m0

2

m,m0 =1

=

M
X

em ⊗ KΦ∗m SCx∗ wm u∗

2

m=1

= K 2 kuk22
= K 2 kuk22
= K 2 kuk22

M
X
m=1
M
X

kΦ∗m SCx∗ wm k22
∗
k(ID ⊗ wm
Cx S ∗ )vec(Φm )k2

m=1
M
X


∗
∗
[vec(Φm )> , . . . , vec(Φm )> ]>
e>
⊗
I
⊗
w
C
S
D
m
m x

2

.

m=1

Therefore, by Lemma A.3,
p
√
k(q)k ≤ C(β)Kkuk2 M Dρx,w log K
holds with probability 1 − K −β .
√
By collecting the estimates with the fact that µ ≤ M , we obtain that with probability 1−K −β ,
the relative perturbation due to (h) is upper bounded by

(54)

√  √
 1
C(β)ρx,w logα (M KL)  √
k(h)k
D
MD
√
≤
µ
M
D
+
K
+
+
K
K 2 kxk22 kuk22
kxk22 kuk2
K 3/2
MK
√




α
0
√
ρx,w
C (β) log (M KL)
K
D
√
√
≤
µ
+√
+ D .
M
ηL
Kσw kxk2
K

Note that the tail estimate in (54) dominates that in (51). It therefore follows that (29) holds
with probability 1 − CK −β . This completes the proof.

SPECTRAL METHODS FOR PASSIVE IMAGING

27

5.4. Proof of Lemma 3.7. The analysis of the noise term En only involves second-order chaos
processes and can be reduced to bounds for suprema of such processes as they are established
in [25, 27].
The first lemma used in the proof bounds the maximum cross-correlation deviation of the noise
terms, which is given by
ρw :=

max

1≤m,m0 ≤M

∗
∗
kS(Cw
Cwm0 − E[Cw
Cwm0 ])S ∗ k.
m
m

See Appendix D.7 for the proof of the lemma.
Lemma 5.9. Suppose that (A2) holds. For any β ∈ N, there is a constant C(β) that depends
only on β such that
√
2
(55)
ρw ≤ C(β)σw
KL logα (M KL)
holds with probability 1 − K −β .
The second proof ingredient is a bound for the average auto-correlation deviation of the noise
terms, which is given by
M
1 X
∗
∗
S(Cw
Cwm − E[Cw
Cwm ])S ∗ .
m
m
M

ρ̄w :=

m=1

The bound is provided by the following lemma; see Appendix D.8 for its proof.
Lemma 5.10. Suppose that (A2) holds. For any β ∈ N, there is a constant C(β) that depends
only on β such that
√
2
(56)
ρ̄w ≤ C(β)σw
M −1/2 KL logα (M KL)
holds with probability 1 − K −β .
For the remainder of the proof of Lemma 3.7, we condition on the event that (55) and (56)
hold.
Define Λ ∈ CM K×M K by
Λ=

M X
M
X

kwm0 k22 em e∗m ⊗ IK .

m=1 m0 =1
m0 6=m

Under the assumption in (A2), we have
2
E[kwm0 k22 ] = σw
L,

∀m0 = 1, . . . , M.

Then it follows that E[Λ] is a scalar multiple of the identity given by
2
E[Λ] = σw
(M − 1)LIM K .

Here we assume that σw is known a priori or can be estimated from the data.
2 (M − 1)LI
We decompose Φ∗ (Yn∗ Yn − σw
M K )Φ into two parts as follows:
(57)

2
Φ∗ (Yn∗ Yn − σw
(M − 1)LIM K )Φ = Φ∗ (Yn∗ Yn − Λ)Φ + Φ∗ (Λ − E[Λ])Φ .
|
{z
} |
{z
}
(r)

(s)

Then we estimate the summands in the right-hand side of (57). The following lemma, which is
proved in Appendix D.9, provides a tail estimate of k(r)k.

28

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Lemma 5.11. Suppose that (A1) holds. For any β ∈ N, there exist a numerical constant α ∈ N
and a constant C(β) that depends only on β such that, conditional on the noise vector w,
√
kΦ∗ (Yn∗ Yn − Λ)Φk ≤ C(β)M KD(ρ̄w + 2ρw ) logα (M D)
holds with probability 1 − K −β .
Next, due to the block diagonal structure, the spectral norm of the second term (s) is upper
bounded by
max

M
X

1≤m≤M

≤



(kwm0 k22 − E[kwm0 k22 ])Φ∗m Φm

m0 =1
m0 6=m
M
X

max

1≤m≤M

kwm0 k22 − E[kwm0 k22 ]



m0 =1
m0 6=m


max kΦ∗m Φm k .

1≤m≤M

2 is a χ2 random variable with LM degrees of freedom. By Lemma A.1,
The first factor divided by σw

max

1≤m≤M

M
X

√
2
LM log(M K)
kwm0 k22 − E[kwm0 k22 ] ≤ C1 (1 + β)σw

m0 =1
m0 6=m

holds with probability 1 − K −β . On the other hand, by (53),
max kΦ∗m Φm k ≤ C2 (1 + β)K log(M K)

1≤m≤M

holds with probability 1 − K −β . Therefore,

√
√
2
k(s)k ≤ C(β)σw
M K L log2 (M K)

with probability 1 − 2K −β .
Finally, (30) follows by collecting the above estimates. This completes the proof.
5.5. Proof of Lemma 3.8. Similarly to the proof of Lemma 3.7, through the triangle inequality,
the left-hand side of (35) is upper bounded by
(58)

kΦ∗ (Yn∗ Yn − Λ)Φuk2 kΦ∗ (Λ − E[Λ])Φuk2
+
.
kuk2
kuk2

For the first term in the right-hand side of (58), we modify Lemma 5.11 as follows; see Appendix D.10 for the proof.
Lemma 5.12. Suppose that (A1) holds. For any β ∈ N, there exist a numerical constant α ∈ N
and a constant C(β) that depends only on β such that, conditional on the noise vector w,
√
√
kΦ∗ (Yn∗ Yn − Λ)Φuk2
≤ C(β) M KD( M ρ̄w + ρw ) logα (M D)
kuk2
holds with probability 1 − K −β .
Combining Lemmas 5.9,5.10, and 5.12 implies √
that the first term in the right-hand side of (58)
is smaller than kΦ∗ (Yn∗ Yn − Λ)Φk by a factor of M .
For the second term in the right-hand side of (58), we use the fact that it is no larger than
∗
kΦ (Λ − E[Λ])Φk. Then we may use the tail estimate derived in the proof of Lemma 3.7.
By collecting the estimates, we obtain that (35) holds with probability 1 − CK −β . This completes the proof.

SPECTRAL METHODS FOR PASSIVE IMAGING

29

6. Conclusion. We studied a passive imaging problem with multiple channels, which is formulated as multichannel blind deconvolution with noise-like source and time-limited impulse responses.
Additionally, motivated by several real world applications, we impose that the FIR coefficients of
impulse responses belong to corresponding low dimensional subspaces. For such a scenario, we proposed a spectral method called subspace-constrained cross-convolution (SCCC) that modifies and
improves upon a classical method developed in the 1990s by overcoming the noise sensitivity. SCCC
provides stable estimates of the impulse responses from finitely many samples and its performance
is backed by theoretical error bounds under generic subspace models. In this scenario, SCCC also
empirically outperforms competing approaches. The fundamental estimates in the analysis of this
paper extend to the sparsity or low-rank cases with minor changes. Corresponding recovery results
on these extended models will be left for follow-up work.
Appendix A. Toolbox.
A.1. Concentration of χ2 Random Variables.
Lemma A.1 (Complexification of [26, Lemma 1]). Let g1 , . . . , gn be independent copies of a
>
standard
Pn complex 2 Gaussian variable. Let a1 , . . . , an be nonnegative and a = [a1 , . . . , an ] . Let
Z = k=1 ak (|gk | − 1). Then, for any t > 0,
√
√
P(Z ≥ 2kak2 t + kak∞ t) ≤ exp(−t),
√
√
P(Z ≤ − 2kak2 t) ≤ exp(−t).
A.2. Hanson-Wright Inequality.
Lemma A.2 (Complexification of [40, Theorem 1.1]). Let A ∈ Cm×n . Let g ∈ Cn be a standard
complex Gaussian vector. For any 0 < ζ < 1, there exists an absolute constant C such that
p
|kAgk22 − E[kAgk22 ]| ≤ C(kA∗ AkF log(2ζ −1 ) ∨ kAk2 log(2ζ −1 ))
holds with probability 1 − ζ.
Lemma A.3 (Complexification of [40, Theorem 2.1]). Let A ∈ Cm×n . Let g ∈ Cn be a standard
complex Gaussian vector. For any 0 < ζ < 1, there exists an absolute constant C such that
p
|kAgk2 − kAkF | ≤ CkAk log(2ζ −1 )
holds with probability 1 − ζ.
A.3. Complexification of Maurey’s Lemma. The following lemma is a direct consequence of
Maurey’s empirical method [5]. Define a block norm on Rmd by
> >
k[q1> , . . . , qm
] k`m
d := max kqk k2 .
∞ (` )
2

1≤k≤m

d
m d
Let `m
∞ (`2 (R)) denote the corresponding Banach space. Similarly, `∞ (`2 (C)) is defined over the
complex scalar field.

Lemma A.4 (Maurey’s empirical method [20, Lemma 3.1]). Let k, m, n ∈ N and T : `k1 (R) →
d
`m
∞ (`2 (R)) be a linear operator. Then
Z ∞q
p
log N (T (B`k (R) ), k · k`m
, t)dt ≤ C (1 + log k)(1 + log m)(1 + log m + log d)kT kop ,
d
∞ (` (R))
0

1

2

where k · kop denotes the operator norm.

30

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Lemma A.4 extends to the complex field case, which is shown in the following corollary.
d
Corollary A.5. Let k, m, n ∈ N and T : `k1 (C) → `m
∞ (`2 (C)) be a linear operator. Then
Z ∞q
p
log N (T (B`k (C) ), k · k`m
(1 + log k)(1 + log m)(1 + log m + log d)kT kop .
d (C)) , t)dt ≤ C
(`
∞
2

1

0

Proof of Corollary A.5. Let ι : C → R2 be a natural map defined by
ι(x) = [Re(x), Im(x)]> ,

∀x ∈ C.

By a slight abuse of notation, we apply ι elementwise to Ck , i.e.
ι([x1 , . . . , xk ]> ) = [Re(x1 ), Im(x1 ), . . . , Re(xk ), Im(xk )]> .
Then we have
ιB`k (C) ⊂

√

1

2B`2k (R)
1

and
ιB`m
= B`m
.
d
2d
∞ (` (C))
∞ (` (R))
2

2

Moreover, ι is bijective. Therefore,
Z ∞q
Z
log N (T (B`k (C) ), k · k`m
, t)dt ≤
d
∞ (` (C))
0

1

2

√
log N (T ( 2B`2k (R) ), k · k`m
, t)dt
2d
∞ (` (R))

∞q

0

1

2

Then the assertion follows from Lemma A.4 with a change of variable in the integral.
A.4. Suprema of Gaussian Processes. We use the following lemma that provides tail estimates
of suprema of first order chaos processes.
Lemma A.6. Let ξ ∈ Cn be a standard Gaussian vector with E[ξξ ∗ ] = In , ∆ ⊂ Cn , and 0 < ζ <
There is an absolute constants C such that
Z ∞p
p
∗
−1
sup |f ξ| ≤ C log(ζ )
log N (∆, k · k2 , t)dt

e−1/2 .

f ∈∆

0

holds with probability 1 − ζ.
Proof. Lemma A.6 is a direct consequence of the moment version of Dudley’s inequality [12, p.
263] and a version of Markov’s inequality [12, Proposition 7.11].
Appendix B. Expectations.
Lemma B.1. Under the assumption in (A1),
∗
E[CΦ
e

m um

2
CΦ
e m um ] = Kkum k2 IL .

Proof of Lemma B.1.
∗
E[CΦ
e

m um

2
∗
CΦ
e m um ] = kum k2 E[CS ∗ g CS ∗ g ],

where g ∈ CK is a standard complex Gaussian vector. Let gk denote the kth entry of g. Since
E[gk gl ] = 0 for all k 6= l, we have
E[CS∗ ∗ g CS ∗ g ]

=

K
X
k=1

This completes the proof.

Ce∗k Cek = KIL .

SPECTRAL METHODS FOR PASSIVE IMAGING

31

Lemma B.2. Under the assumption in (A1),
∗
E[CΦ
e

m um

e m ] = Ke1 u∗ .
Φ
m

e m um and Φ
e m q are independent for any q ∈ CK
Proof of Lemma B.2. By (A1), it follows that Φ
∗
such that q um = 0. Therefore,
∗
e
E[CΦ
e m um Φm ]

∗
e
um u∗m i E[CΦ
e m um Φm um ] ∗
e
=
um
Φm
m um
kum k22
kum k22
C e m um e1 ]
E[C ∗e
Φm um Φ
u∗m = Ke1 u∗m ,
=
kum k22

h
∗
= E CΦ
e

where the last step follows from Lemma B.1.
Lemma B.3. Under the assumption in (A1),
(
K 2 kxk22 kum0 k22 ID
∗
e
e ∗ C∗
E[Φ
e 0 u 0 Φm ] =
m Φ
e 0 u 0 Cx Cx CΦ
m
m
m
m
K 2 kxk22 (kum0 k22 ID + um0 u∗m0 )

m 6= m0 ,
m = m0 .

Proof of Lemma B.3. Suppose that m 6= m0 . Then, by the independence of Φm and Φm0 , it
follows that
i
h
i
h
∗
∗
∗
∗
e
e
e ∗ C∗
e
C
C
C
C
C
C
Φ
]
=
E
Φ
E[Φ
Φ
E
C
m
Φm
e 0u 0
e 0u 0 m
x x Φ
x x Φ
m Φ
m Φ m0
e u
e u
Φ
m0

m0

m

m0

m

=
=
=

m0

m

m

e∗ Φ
e
Kkum0 k22 EΦm [Φ
m m]
Kkum0 k22 EΦm [Φ∗m SS ∗ Φm ]
Kkum0 k22 tr(SS ∗ )ID = K 2 kum0 k22 ID ,

where the second identity holds by Lemma B.1. This proves the first case.
Next we assume that m0 = m. For notational simplicity, let P ∈ CD×D denote the orthogonal
projection onto the span of um , i.e.
um u∗m
P =
.
kum k22
Then by (A1) it follows that Φm um and Φm (ID − P ) are independent. Therefore
e ∗ C∗
E[Φ
m Φ
e

m um

e
Cx∗ Cx CΦ
e m um Φm ]

e ∗ C∗
= E[(P + ID − P )Φ
m Φ
e

m um

e
Cx∗ Cx CΦ
e m um Φm (P + ID − P )]

∗
e ∗ C∗
e
= E[P Φ
e m um Φm P ]
m Φ
e m u m Cx Cx CΦ
{z
}
|
(?)

+ E[(ID −
|

∗
e ∗ C∗
e
P )Φ
e m um Φm (ID
m Φ
e m u m Cx Cx CΦ

{z

(??)

− P )] .
}

The first term (?) is rewritten as
E[g ∗ SCS∗ ∗ g Cx∗ Cx CS ∗ g S ∗ gum u∗m ],

32

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

where g ∈ CK is a standard complex Gaussian vector. On the other hand, we have
E[g

∗

SCS∗ ∗ g Cx∗ Cx CS ∗ g S ∗ g]

K
X

=
=

E[gj1 gj2 gj3 gj4 ]e∗j1 Ce∗j Cx∗ Cx Cej3 ej4
2

j1 ,j2 ,j3 ,j4 =1
2Kkxk22 +

2K(K − 1)kxk22 = 2K 2 kxk22 .

Therefore,
(?) = 2K 2 kxk22 um u∗m .
By the independence of Φm um and Φm (ID −P ) together with the commutativity of convolution,
the second term (??) is computed as
e ∗ E[C ∗
(??) = E[(ID − P )Φ
m
e
Φ

m um

∗
e
CΦ
e m um ]Cx Cx Φm (ID − P )]

e ∗ C ∗ Cx Φ
e m (ID − P )]
= Kkum k22 E[(ID − P )Φ
m x
= Kkum k22 tr(SCx∗ Cx S ∗ )(ID − P )
= K 2 kxk22 (kum k22 ID − um u∗m ).
Collecting these results proves the second case and the proof is done.
Appendix C. Proof of Lemmas in Section 3.
C.1. Proof of Lemma 3.9. By the homogeneity of ρx,w , we may assume that σw = 1. Let
e∗ q,
az,q = CS∗e∗ z Cx S
e is defined in (18). Then
where S
e ∗ Cwm S
e∗ k =
kSC
x

sup
z,q∈B23K

|a∗z,q wm |.

We use Lemma A.6 to get a tail estimate.
Since
kaz,q − az0 ,q0 k2 ≤ kaz,q − az,q0 k2 + kaz,q0 − az0 ,q0 k2
e∗ kkq − q 0 k2 + kC e∗ 0 kkC ∗ S
e∗ kkz − z 0 k2
≤ kCS∗e∗ z kkCx S
x
S q
p
≤ Kρx (kq − q 0 k2 + kz − z 0 k2 ),
it follows that the γ2 functional of the set {az,q | z, q ∈ B23K } is upper bounded by
Z ∞q
p
√
C1 Kρx
log N (B23K , k · k2 , t)dt ≤ C2 K ρx ,
0

where we used a standard volume argument. The assertion follows from Lemma A.6 and a union
bound argument.
C.2. Proof of Lemma 3.10. By the homogeneity, we may assume that σx = σw = 1. Then
e ∗ Cx S
e∗ k = sup |x∗ C ∗ C e∗ x|
kSC
x
e∗ z S z
S
z∈B23K

SPECTRAL METHODS FOR PASSIVE IMAGING

33

and
e∗ k =
e ∗ Cwm S
kSC
x

sup
z,q∈B23K

|x∗ CS∗e∗ q CSe∗ z wm |

are suprema of second order chaos processes. We estimate their tail decay using Theorems 5.1 and
5.2. (For more details, see the proof of Lemma 5.9.) By Theorem 5.1 and a union bound argument,
it follows that
√
ρx ≤ Cβ σx2 KL log5 (KL)
holds with probability 1 − K −β . Similarly, by Theorem 5.2,
√
ρx,w ≤ Cβ σw σx KL log5 (M KL)
holds with probability 1 − K −β .
Moreover, by Lemma A.1, we also have that
kxk22 ≥ σx2 (L −

p
2Lβ log K)

holds with probability 1 − K −β . The assertions follow by assembling the above estimates.
Appendix D. Proof of Lemmas in Section 5.
D.1. Proof of Lemma 5.3. First note that kΨ∗ AΨ − E[Ψ∗ AΨ]k is written as a variational
form given by
max |q ∗ (Ψ∗ AΨ − E[Ψ∗ AΨ])q 0 |.

(59)

q,q 0 ∈B2D

For all q ∈ CD , we have
Ψq = (q > ⊗ IK )vec(Ψ).
Let ψ = vec(Ψ). Then
q ∗ Ψ∗ AΨq 0 = ψ ∗ (q > ⊗ IK )∗ A((q 0 )> ⊗ IK )ψ = ψ ∗ (q(q 0 )> ⊗ A)ψ.
Then (59) becomes the supremum of the second order chaos process
ψ ∗ (q(q 0 )> ⊗ A)ψ − E[ψ ∗ (q(q 0 )> ⊗ A)ψ]
over q ∈ B2D . We obtain its tail estimate by applying Theorem 5.2 with
∆1 = {q > ⊗ IK | q ∈ B2D }
and
∆2 = {(q 0 )> ⊗ A | q 0 ∈ B2D }.
By direct calculation, we obtain
dS (∆1 ) ≤ 1,

dF (∆1 ) ≤

dS (∆2 ) ≤ kAk,

dF (∆2 ) ≤

and

√
K,
√
KkAk.

34

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Moreover
∞q

Z
γ2 (∆1 , k · k) ≤ C1
0

N (B2D , k

Similarly, we also have

Z

1

r

· k2 , t)dt ≤ C1
0


√
2
2D log 1 +
dt ≤ C2 D.
t

√
γ2 (∆2 , k · k) ≤ C2 DkAk.

The assertion follows from Theorem 5.2 with
s
γ2 (∆2 , k · k)dF (∆2 )
(60)
a=
.
γ2 (∆1 , k · k)dF (∆1 )
D.2. Proof of Lemma 5.4. Note that kΥs,m − E[Υs,m ]k is written as a variational form given
by
sup |z ∗ (Υs,m − E[Υs,m ])z|.
z∈B2K

By the commutativity of convolution, we have
∗

z Υs,m z =

M
X

z ∗ SCS∗ ∗ Φm0 um0 Cx∗ Cx CS ∗ Φm0 um0 S ∗ z

m0 =1
m0 6=m

=

M
X

u∗m0 Φ∗m0 SCS∗ ∗ z Cx∗ Cx CS ∗ z S ∗ Φm0 um0

m0 =1
m0 6=m

=

M
X

vec(Φm0 )∗ (um0 ⊗ IK )SCS∗ ∗ z Cx∗ Cx CS ∗ z S ∗ (u>
m0 ⊗ IK )vec(Φm0 )

m0 =1
m0 6=m

=

M
X

∗
∗
∗
vec(Φm0 )∗ (um0 u>
m0 ⊗ SCS ∗ z Cx Cx CS ∗ z S )vec(Φm0 ),

m0 =1
m0 6=m

where the third identity follows from
Φm0 um0 = (u>
m0 ⊗ IK )vec(Φm0 ).

(61)
Let

Q(z) =

M
X

∗
∗
em0 e∗m0 ⊗ u>
m 0 ⊗ Cx CS ∗ z S ,

m0 =1
m0 6=m

and
φ = [vec(Φ1 )> , . . . , vec(ΦM )> ]> .
Then φ follows the distribution CN (0M KD,1 , IM KD ) and
z ∗ Υs,m z = φ∗ Q(z)∗ Q(z)φ.
Therefore,
sup |z ∗ (Υs,m − E[Υs,m ])z| = sup |φ∗ Q(z)Q(z)φ|.
z∈B2K

z∈B2K

SPECTRAL METHODS FOR PASSIVE IMAGING

35

We get a tail bound of the supremum of the second order chaos process φ∗ Q(z)∗ Q(z)φ −
E[φ∗ Q(z)∗ Q(z)φ] by applying Theorem 5.1 with M = Q(z) and ∆ = {Q(z)|z ∈ B2K }.
Recall that in Section 5.2 we defined S̆ by


0K−1,L−K+1
IK−1
S̆ =
.
IK
0K,L−K
Then the radius of ∆ with respect to the Frobenius norm is upper bounded by
dF (∆) ≤ sup kuk2 kCx S̆ ∗ kkSCS ∗ z kF
z∈B2K

p
√
≤ sup kuk2 kS̆Cx∗ Cx S̆ ∗ k1/2 Kkzk2 ≤ kuk2 ρx K,
z∈B2K

where the first inequality follows from the identity CS∗ ∗ z S ∗ = S̆ ∗ S̆CS∗ ∗ z S ∗ .
Let k · k2,∞ be defined in (39). Then the radius of ∆ with respect to the spectral norm is upper
bounded by
dS (∆) ≤ sup kuk2,∞ kCx S̆ ∗ kkSCS ∗ z k
z∈B2K

√ √
≤ sup kuk2,∞ ρx LkF S ∗ zk∞

(62)

z∈B2K

p
√
≤ sup kuk2,∞ ρx kzk1 ≤ kuk2,∞ ρx K,
z∈B2K

√
where the second inequality follows from the identity CS ∗ z = LF ∗ diag(F S ∗ z)F and the third
L
√1
inequality follow from kF : `L
1 → `∞ k = L .
By (62) and Dudley’s inequality, the γ2 functional of ∆ is upper bounded by
Z ∞q
√ √
γ2 (∆, k · k) ≤ C1 kuk2,∞ ρx L
log N (F S ∗ B2K , k · k∞ , t)dt
0
√ √ p
≤ C2 kuk2,∞ ρx K log K log3/2 L,
where the last step follows from Corollary A.5. The assertion follows by applying these estimates
to Theorem 5.1 together with a union bound argument over m = 1, . . . , M .
D.3. Proof of Lemma 5.5. Note that kZm − E[Zm ]k is written as a variational form given by
(63)

sup

|z ∗ Zm q − E[z ∗ Zm q]|.

z∈B22L−1 ,q∈B2D

Let φm = vec(Φm ). Then
z ∗ Zm q = z ∗ S̆CS∗ ∗ Φm um S ∗ Φm q = φ∗m (um q > ⊗ SCS̆∗ ∗ z S ∗ )φm .
Then (63) becomes the supremum of the second order chaos process
∗
∗ ∗ >
∗
>
∗
∗ ∗ >
φ∗m (u>
m ⊗ SCS̆ ∗ z S ) (q ⊗ IK )φm − E[φm (um ⊗ SCS̆ ∗ z S ) (q ⊗ IK )φm ]

over z ∈ B22K−1 and q ∈ B2D . We obtain its tail estimate by applying Theorem 5.2 with
2K−1
∗
∗
∆1 = {u>
}
m ⊗ SCS̆ ∗ z S | z ∈ B2

36

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

and
∆2 = {q > ⊗ IK | q ∈ B2D }.
The radii of ∆1 and ∆2 are upper bounded by
√
dS (∆1 ) ≤ dF (∆1 ) ≤ kuk2,∞ K
and
dS (∆2 ) ≤ 1,
Moreover, since

√

dF (∆2 ) ≤

K.

√
∗
∗
∗
ku>
m ⊗ SCS̆ ∗ z S k ≤ kuk2,∞ LkF S̆ zk∞ ,

we have
√ Z
γ2 (∆1 , k · k) ≤ C1 kuk2,∞ L

∞q

N (F S̆ ∗ B22K−1 , k · k∞ , t)dt
0
Z ∞q
√
≤ C2 kuk2,∞ LK
N (F S̆ ∗ B12K−1 , k · k∞ , t)dt
√ p 0
≤ C3 kuk2,∞ K log(2K − 1) log3/2 L,

where the last step follows from Corollary A.5.
On the other hand, we also have
Z ∞q
Z
D
γ2 (∆2 , k · k) ≤ C1
N (B2 , k · k2 , t)dt ≤ C1
0

1

0

r


√
2
2D log 1 +
dt ≤ C4 D.
t

By applying the estimates to Theorem 5.2 with a given in (60), we obtain that the supremum
is upper bounded by
p
√ p
√
C(kuk2,∞ K 3/2 D log K log3/2 L + kuk2,∞ K log K log3/2 L + kuk2,∞ KD) log(log M + β log K)
with probability 1 − M −1 K −β .
The assertion follows by applying a union bound argument over m = 1, . . . , M .
P
∗
D.4. Proof of Lemma 5.6. Note that the spectral norm of M
m=1 em ⊗ (Zm − E[Zm ]) admits
a variational form given by
(64)

sup
z∈B22L−1 ,v∈B2M D

M
X

z ∗ Zm vm − E[z ∗ Zm vm ] ,

m=1

> ]> .
where vm ∈ CD for m = 1, . . . , M and v = [v1> , . . . , vM
>
>
Let φm = vec(Φm ) for m = 1, . . . , M and φ = [φ1 , . . . , φ>
M ] . Then as before
M
X

∗

z Zm vm =

m=1

=

M
X
m=1
M
X

z ∗ S̆CS∗ ∗ Φm um S ∗ Φm vm
>
φ∗m (um vm
⊗ SCS̆∗ ∗ z S ∗ )φm

m=1

= φ∗

M
X
m=1

∗
∗
em e∗m ⊗ u>
m ⊗ SCS̆ ∗ z S

M
∗  X
m=1


>
em e∗m ⊗ vm
⊗ IK φ.

SPECTRAL METHODS FOR PASSIVE IMAGING

37

Then (64) becomes the supremum of a second order chaos process. We obtain its tail estimate
by applying Theorem 5.2 with
∆1 =

M
nX

2K−1
∗
∗
em e∗m ⊗ u>
m ⊗ SCS̆ ∗ z S | z ∈ B2

o

m=1

and
∆2 =

M
nX

o
>
em e∗m ⊗ vm
⊗ IK | v ∈ B2M D .

m=1

Let k · k2,∞ be defined in (39). Then the radii of ∆1 and ∆2 are upper bounded by
√
√
√
dS (∆1 ) ≤ kuk2,∞ K, dF (∆1 ) ≤ kuk2 K ≤ kuk2,∞ M K,
and
dS (∆2 ) ≤ 1,

dF (∆2 ) ≤

√
K.

Moreover, since
M
X

√
∗
∗
≤ kuk2,∞ LkF S̆ ∗ zk∞ ,
em e∗m ⊗ u>
m ⊗ SCS̆ ∗ z S

m=1

we have
√ Z
γ2 (∆1 , k · k) ≤ C1 kuk2,∞ L

∞q

N (F S̆ ∗ B22K−1 , k · k∞ , t)dt
0
Z ∞q
√
≤ C2 kuk2,∞ LK
N (F S̆ ∗ B12K−1 , k · k∞ , t)dt
0
√ p
≤ C3 kuk2,∞ K log(2K − 1) log3/2 L,

where the last step follows from Corollary A.5.
On the other hand, since
M
X

>
em e∗m ⊗ vm
⊗ IK ≤ kvk2,∞ ,

m=1

we also have
Z
γ2 (∆2 , k · k) ≤ C1
0

√ p
N (B2M D , k · k2,∞ , t)dt ≤ C4 D log D log(M D).

∞q

By applying these estimates to Theorem 5.2 with a given in (60), we obtain that the supremum
is upper bounded by
√
C 0 (β) logα (M KL)kuk2,∞ (M 1/4 K 3/4 D1/4 + K + M KD)
with probability 1 − K −β . Finally, by the inequality of arithmetic and geometric means, we have
√
K + M KD
1/4 3/4 1/4
M K D
≤
.
2
This completes the proof.

38

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

D.5. Proof of Lemma 5.7. Similar to the proof of Lemma 5.4, we rewrite kΥc,m k as a variational form given by
sup |z1∗ Υc,m z2 |.
z1 ,z2 ∈B2K

By the commutativity of convolution and (61), z1∗ Υc,m z2 is written as follows
M
X

z1∗ SCS∗ ∗ Φm0 um0 Cx∗ Cwm0 S ∗ z2

m0 =1
m0 6=m

M
X

=

vec(Φm0 )∗ (um0 ⊗ SCS ∗ z1 Cx∗ Cwm0 S ∗ z2 ),

m0 =1
m0 6=m

which, conditional on w, is a centered Gaussian process indexed by (z1 , z2 ) ∈ B2K × B2K . We
compute a tail estimate of the supremum by applying Lemma A.6 with
∆ = {f (z1 , z2 ) | z1 , z2 ∈ B2K },
where
f (z1 , z2 ) =

M
X

em0 ⊗ um0 ⊗ SCx∗ Cwm0 CS ∗ z1 S ∗ z2 .

m0 =1
m0 6=m

Since kCS ∗ (zj −zj0 ) k =

√

LkF S ∗ (zj − zj0 )k∞ for j = 1, 2 and hence

kf (z1 , z2 ) − f (z10 , z20 )k2 ≤ kf (z1 , z2 ) − f (z10 , z2 )k2 + kf (z10 , z2 ) − f (z10 , z20 )k2
√
≤ kuk2 ρx,w L(kF S ∗ (z1 − z10 )k∞ + kF S ∗ (z2 − z20 )k∞ ),
it follows that
Z ∞p
√ Z
log N (∆, k · k2 , t)dt ≤ C1 kuk2 ρx,w L
0

∞q

log N (F S ∗ B2K , k · k∞ , t)dt
0
Z ∞q
√
≤ C1 kuk2 ρx,w KL
log N (F S ∗ B1K , k · k∞ , t)dt
0
√ p
≤ C2 kuk2 ρx,w K log K log3/2 L,

where the last step follows from Corollary A.5 together with the observation that kF k`L →`L∞ ≤
1

L−1/2 . Then the assertion follows from Lemma A.6 with a union bound argument over m =
1, . . . , M .
>
D.6. Proof of Lemma 5.8. Let φm = vec(Φm ) for m = 1, . . . , M and φ = [φ>
, . . . , φ>
1P
M] .
M
> ] ∈ B M D . Then the spectral norm of
∗
Let v1 , . . . , vM ∈ CD and v = [v1> , . . . , vM
2
m=1 em ⊗
∗
∗
S̆Cwm Cx S Φm is rewritten as
M
X

∗
e∗m ⊗ S̆Cw
Cx S ∗ Φm =
m

m=1

=

=

sup

sup

z∈B22K−1

v∈B2M D

sup

sup

z∈B22K−1

v∈B2M D

sup
z∈B22K−1

sup

M
X

∗
z ∗ S̆Cw
Cx S ∗ Φm vm
m

m=1
M
X

>
∗
(vm
⊗ z ∗ S̆Cw
Cx S ∗ )φm
m

m=1
M
X

v∈B2M D m=1

>
∗
e∗m ⊗ (vm
⊗ z ∗ S̆Cw
Cx S ∗ )φ .
m

SPECTRAL METHODS FOR PASSIVE IMAGING

Let
f (z, v) =

M
X

39

em ⊗ (vm ⊗ SCx∗ Cwm S̆ ∗ z).

m=1

Then we obtain
M
X

∗
e∗m ⊗ S̆Cw
Cx S ∗ Φm = sup
m

v∈B2M D

m=1

|f (z, v)∗ φ|.

sup
z∈B22K−1

Note that f (z, v)∗ φ, conditioned on w, is a centered Gaussian process. We compute a tail estimate
of this supremum by applying Lemma A.6 with
∆ = {f (z, v) | z ∈ B22K−1 , v ∈ B2M D }.
Then we need to compute the entropy integral for ∆. Recall
e ∗ Cw S
e∗ k ≥ kS̆C ∗ Cx S ∗ k,
ρx,w = max kSC
x
wm
1
1≤m≤M

∀m = 1, . . . , M.

By the triangle inequality, we obtain
kf (z, v) − f (z 0 , v 0 )k2 ≤ kf (z, v) − f (z, v 0 )k2 + kf (z, v 0 ) − f (z 0 , v 0 )k2
≤ ρx,w (kzk2 kv − v 0 k2 + kz − z 0 k2 kv 0 k2 )
≤ ρx,w (kv − v 0 k2 + kz − z 0 k2 ).
The integral of the log-entropy number is computed as
sup

sup

v∈B2M D

z∈B22K−1

Z
≤ C1

|f (z, v)∗ φ|

∞p

log N (∆, k · k2 , t)dt
Z ∞q

∞q
M
D
log N (B2 , k · k2 , t)dt +
log N (B22K−1 , k · k2 , t)dt
≤ C1 ρx,w
0
√0
√
≤ C2 ρx,w ( M D + K),
0

Z

where the last step follows from a standard volume argument. Then the assertion follows from
Lemma A.6.
D.7. Proof of Lemma 5.9. By the homogeneity of ρw , we may assume that σw = 1. We first
consider the case that m0 = m. Let z ∈ CK . Then by the commutativity of convolution, we have
Cwm S ∗ z = wm ~ S ∗ z = CS ∗ z wm .
∗ C
∗
∗
Therefore, the spectral norm of S(Cw
wm − E[Cwm Cwm ])S is rewritten as
m
∗
∗
∗
∗
kS(Cw
Cwm − E[Cw
Cwm ])S ∗ k = sup |z ∗ S(Cw
Cwm − E[Cw
Cwm ])S ∗ z|
m
m
m
m
z∈B2K

∗
∗
CS∗ ∗ z CS ∗ z wm ]|,
= sup |wm
CS∗ ∗ z CS ∗ z wm − E[wm
z∈B2K

where the last term is the supremum of a second order chaos. We use Theorem 5.1 to get its tail
estimate. Define
∆z = {CS ∗ z | z ∈ B2K }.

40

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Then the radii of ∆z with respect to the spectral and Frobenius norms are given by
√
√
dS (∆z ) = sup kCS ∗ z k = sup LkF S ∗ zk∞ = sup kzk1 = K
z∈B2K

z∈B2K

dF (∆z ) = sup kCS ∗ z kF = sup
z∈B2K

√

z∈B2K

√
Lkzk2 =

L.

z∈B2K

Moreover, the γ2 functional of ∆z is bounded by
√ Z ∞q
log N (F S ∗ B2K , k · k∞ , t)dt
γ2 (∆z , k · k) ≤ C1 L
0
Z ∞q
√
log N (F S ∗ B1K , k · k∞ , t)dt
≤ C1 KL
0
√ p
≤ C1 K log K log3/2 L,
where the last step follows from Lemma A.4. By Theorem 5.1,
∗
∗
2
kS(Cw
Cwm − E[Cw
Cwm ])S ∗ k ≤ Cσw
K(1 + log K) log3 L(1 + 2 log M + β log K)
m
m

holds with probability 1 − M −2 K −β .
∗ C
Next we consider the case where m0 =
6 m. In this case, we have E[Cw
wm0 ] = 0L,L . Similarly
m
∗
∗
to the previous case, the spectral norm of SCwm Cwm0 S is rewritten as
∗
∗
∗
∗
em,m
em,m0 |,
kSCw
Cwm0 S ∗ k = sup |wm
CS∗ ∗ q CS ∗ z wm0 | = sup |w
0 Lq Rz w
m
z,q∈B2K

where



em,m0
w


wm
=
,
wm 0

z,q∈B2K




0L,L 0L,L
Lq =
,
CS ∗ q 0L,L




0L,L CS ∗ z
Rz =
.
0L,L 0L,L

Define
∆Lq = {Lq | q ∈ B2K }
and
∆Rz = {Rz | z ∈ B2K }.
Then, the radii and γ2 functional of ∆Lq and ∆Rz are identical to those of ∆z . Therefore, by
Theorem 5.2,
∗
∗
2
kS(Cw
Cwm0 − E[Cw
Cwm0 ])S ∗ k ≤ Cσw
K(1 + log K) log3 L(1 + 2 log M + β log K)
m
m

holds with probability 1 − M −2 K −β . The assertion follows by applying a union bound argument.
D.8. Proof of Lemma 5.10. Without loss of generality, we assume that σw = 1. Similarly
to proof of Lemma 5.9 in Appendix D.7, we can rewrite M ρ̄w as the supremum of a second order
chaos process as follow:
M ρ̄w = sup
z∈B2K

M
X

∗
∗
wm
CS∗ ∗ z CS ∗ z wm − E[wm
CS∗ ∗ z CS ∗ z wm ]

m=1

= sup w∗ (IM ⊗ CS∗ ∗ z CS ∗ z )w − E[w∗ (IM ⊗ CS∗ ∗ z CS ∗ z )w] ,
z∈B2K

> ]> is a standard complex Gaussian vector of length M L.
where w = [w1> , . . . , wM

SPECTRAL METHODS FOR PASSIVE IMAGING

41

Define
e z = {IM ⊗ CS ∗ z | z ∈ B K }.
∆
2
e z with respect to the spectral and Frobenius norms are upper bounded respecThen the radii of ∆
tively by
√
√
e z ) = sup LkF S ∗ zk∞ ≤ K
dS (∆
z∈B2K

e z ) = sup
dF (∆

√

M kCS ∗ z kF ≤

√
M L.

z∈B2K

e z is upper bounded by
Moreover, since kIM ⊗ CS ∗ z k = kCS ∗ z k, the γ2 functional of ∆
√ p
e z , k · k) = γ2 (∆z , k · k) ≤ C1 K log K log3/2 L,
γ2 (∆
where the last step has been shown in Appendix D.7. By Theorem 5.1,
√
2
ρ̄w ≤ Cβσw
M −1/2 KL(1 + log K) log3 L log K
holds with probability 1 − K −β .
D.9. Proof of Lemma 5.11. First we rewrite Φv as
Φv = Wv φ,
where
φ = [vec(Φ1 )> , . . . , vec(Φ1 )> ]>
and
Wv =

M
X

>
em e∗m ⊗ vm
⊗ IK .

m=1

Then it follows that
v ∗ Φ∗ (Yn∗ Yn − Λ)Φv = φ∗ Wv∗ (Yn∗ Yn − Λ)Wv φ,
where the latter, conditional on w, is a quadratic Gaussian form. Furthermore, by direct calculation,
we have
Eφ [φ∗ Wv∗ (Yn∗ Yn − Λ)Wv φ] = Eφ [Φ∗ (Yn∗ Yn − Λ)Φ] = tr(Yn∗ Yn − Λ)IM D = 0M D,M D .
Then kΦ∗ (Yn∗ Yn − Λ)Φk is written as
sup |φ∗ Wv∗ (Yn∗ Yn − Λ)Wv φ − Eφ [φ∗ Wv∗ (Yn∗ Yn − Λ)Wv φ]|,

v∈B2M D

which is the supremum of a second order Gaussian chaos process. We compute its tail estimate by
applying Theorem 5.2 with
∆1 = {Wv | v ∈ B2M D }
and
∆2 = {(Yn∗ Yn − Λ)Wv | v ∈ B2M D }.

42

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Let k · kp,q be defined in (39). Then the radii of ∆1 with respect to the Frobenius and spectral
norms are upper bounded respectively by
dS (∆1 ) ≤ sup kvk2,∞ ≤ 1
v∈B2M D

dF (∆1 ) ≤ sup

√

Kkvk2 ≤

√
K.

v∈B2M D

By Lemma E.1, the γ2 functional is bounded by
Z ∞q
√ p
γ2 (∆1 , k · k) ≤ C1
log N (B2M D , k · k2,∞ , t)dt ≤ C2 D log D log(M D).
0

We repeat the calculation for ∆2 . Note that (Yn∗ Yn − Λ)Wv is expressed as
(Yn∗ Yn

− Λ)Wv =

M
X

em e∗m

⊗

>
vm

⊗

M X
M
X
m=1

∗
S(Cwm0 Cw
− kwm0 k22 IL )S ∗
m0

m0 =1
m0 6=m

m=1

−

M
X

>
∗
∗
em e∗m0 ⊗ vm
0 ⊗ SCwm Cw 0 S .
m

m0 =1
m0 6=m

Noting that all the summands in this decomposition are K × K matrices and orthogonal with
respect to the Frobenius inner product, we obtain that the Frobenius norm of (Yn∗ Yn − Λ)Wv is
upper bounded by
k(Yn∗ Yn

−

Λ)Wv k2F

≤

M
X
m=1

+

M
X

kvm k22 K

M X
M
X

2

m0 =1
m0 6=m
∗
kvm0 k22 KkSCwm Cw
S ∗ k2
m0

m=1 m0 =1
m0 6=m
≤ 2kvk22 K(M 2 ρ̄2w

which implies that

∗
S(Cwm0 Cw
− kwm0 k22 IL )S ∗
m0

+ ρ2w ) + kvk22 KM ρ2w ,

√
√
k(Yn∗ Yn − Λ)Wv kF ≤ Ckvk2 K(M ρ̄w + M ρw ).

On the other hand, it follows from the block Gershgorin disk theorem [11] that
k(Yn∗ Yn − Λ)Wv k

≤ max kvm k2
1≤m≤M

M
X

∗
S(Cwm0 Cw
− kwm0 k22 IL )S ∗ +
m0

m0 =1
m0 6=m

M
X

∗
kvm0 k2 kSCwm Cw
S∗k
m0



m0 =1
m0 6=m

≤ kvk2,∞ (M ρ̄w + ρw ) + kvk2,1 ρw ≤ kvk2,∞ M (ρ̄w + 2ρw ).
Therefore, the radii of ∆2 with respect to the Frobenius and spectral norms are upper bounded
respectively by
dS (∆2 ) ≤ M (ρ̄w + 2ρw ),
√
√
dF (∆2 ) ≤ C K(M ρ̄w + M ρw ).

SPECTRAL METHODS FOR PASSIVE IMAGING

43

Moreover, by Lemma E.1, the γ2 functional is bounded by
Z ∞q
log N (B2M D , k · k2,∞ , t)dt
γ2 (∆2 , k · k) ≤ C1 M (ρ̄w + 2ρw )
0
p
√
≤ C2 M D(ρ̄w + 2ρw ) log D log(M D).
Consequently, Theorem 5.2 yields that
√
kΦ∗ (Yn∗ Yn − Λ)Φk ≤ C(β)M KD(ρ̄w + 2ρw ) log D log2 (M D)
holds with probability at least 1 − K −β .
D.10. Proof of Lemma 5.12. We modify the proof of Lemma 5.11 in Appendix D.9 as follows.
First note that
|φ∗ Wv∗ Q∗w Qw Wu φ|
kΦ∗ (Yn∗ Yn − Λ)Φuk2
= sup
.
kuk2
kuk2
v∈B M D
2

˘ R by the following singleton set
We only need to replace ∆
˘ R,u = {(Yn∗ Yn − Λ)Wu }.
∆
˘ R,u and ∆
˘ R are the same for both the Frobenius and spectral norms. HowIndeed, the radii of ∆
˘
ever, the γ2 functional of ∆R,u is 0. The assertion follows by applying the modified estimates to
Theorem 5.2.
Appendix E. Entropy Estimates.
Lemma E.1.
Z
(65)
0

√ p
log N (B2M D , k · k2,∞ , t)dt ≤ C D log D log(M D).

∞q

Proof of Lemma E.1. Let us recall that the (2, ∞)-block norm of v ∈ CM D is defined by
kvk2,∞ = max kvm k2 ,
m∈[M ]

> ].
where vk ∈ CD for k = 1, . . . , M denotes the blocks of v such that v = [v1> , . . . , vM
Since
kvk2,∞ ≤ kvk2 ,

the interval in the integral in (65) can be restricted to the unit interval [0, 1).
Indeed, the (2, ∞)-block norm of v is rewritten as
kvk2,∞ = max max hvm , ζi.
m∈[M ] ζ∈B2D

To compute an estimate of the entropy integral in (65), we adopt the strategy [20] that estimates
a unit ball using a polytope. The original strategy [20] was developed for the RIP analysis for lowrank tensors and applied to the tensor nuclear norm. The same strategy applies to the block norm
in this section too.
Lemma E.2. There exist ζ1 , . . . , ζN ∈ SD−1 such that
B2D ⊂ 2 conv{ζ1 , . . . , ζN }
and log N ≤ (2D + 1) log(4D + 3).

44

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

Proof of Lemma E.2. Let {ζ1 , . . . , ζN } ⊂ SD−1 be an -net of SD−1 . Then by the standard
volume argument, we have N ≤ (1 + 2/)2D . Furthermore, it follows that
SD−1 ⊂

1
· absconv{ζ1 , . . . , ζN }.
1−

bk , αk )k∈N ⊂ {ζ1 , . . . , ζN } × C as follows. Let
Indeed, for any w ∈ SD−1 , we construct a sequence (w
b1 be the closest vector to w in {ζ1 , . . . , ζN }. If w = α1 w
b1 , then αk = 0 and w
bk = 0
α1 = 1 and w
b1 k2 and w
b2 be the closest vector to α2−1 (w − α1 w
b1 )
for all k ≥ 2. Otherwise, let α2 = kw − α1 w
bk = 0 for all k ≥ 3. Otherwise, let
b2 , then αk = 0 and w
b1 + α2 w
in {ζ1 , . . . , ζN }. If w = α1 w
b2 ) in {ζ1 , . . . , ζN }.
b 2 − α1 w
b2 k2 and w
b3 be the closest vector to α3−1 (w − α1 w
b1 − α2 w
α3 = kw − α1 w
By continuing in this way, we have
X
bk ,
w=
αk w
k∈N

bk k2 = 1 for all k ∈ N. Therefore,
where |αk | ≤ k−1 and kw
X

|αk | ≤

k∈N

1
1−

and the assertion follows. By including ±ζk instead of ζk , we can replace the absolute convex hull
by convex hull and the cardinality increases only by factor 2. Choosing  = 1/(2D + 1) completes
the proof.
By Lemma E.2, we approximate the (2, ∞)-block norm of v as a polytope norm as follows:
kvk2,∞ = max max hvm , ζi
ζ∈SD−1 m∈[M ]

≤2
(66)

max

max hvm , ζi

ζ∈conv{ζ1 ,...,ζN } m∈[M ]

= 2 max max |hvm , ζn i|
n∈[N ] m∈[M ]

=: 2 ||| v ||| ,
where log N ≤ (D + 1) log(2D + 3).
Define
n
> >
Sζn = v = [v1> , . . . , vM
]

o
max |hvm , ζn i| ≤ 1 .

m∈[M ]

Then its polar set is given by
Sζ◦n = {z ⊗ ζn |z ∈ B1M } = conv{ed ⊗ ζn |d ∈ [D]}.
T
Note that the unit ball with respect to ||| · ||| is given as n∈[N ] Sζn . To compute the unit ball
with respect to the dual norm, we will use a well known polar duality result in the following lemma.
Note that Sζn is not bounded. As the lemma is typically stated for bounded sets, we provide the
proof for completeness, verifying that boundedness is not a crucial assumption.
Lemma E.3. Let A and B be convex sets. The polar set of the intersection of A and B is given
by
(A ∩ B)◦ = conv(A◦ ∪ B ◦ ).

SPECTRAL METHODS FOR PASSIVE IMAGING

45

Proof of Lemma E.3. We first show
(conv(A◦ ∪ B ◦ ))◦ ⊂ A ∩ B.

(67)

Suppose that x 6∈ A ∩ B. Since A and B are convex, we have A = A◦◦ and B = B ◦◦ . Without
loss of generality, we may assume that x 6∈ A◦◦ . Then there exists w ∈ A◦ such that
|hx, wi| > 1.

(68)

Since A◦ ⊂ A◦ ∪ B ◦ ⊂ conv(A◦ ∪ B ◦ ), w also satisfies
w ∈ conv(A◦ ∪ B ◦ ).

(69)

The existence of w satisfying both (68) and (69) implies that
x 6∈ (conv(A◦ ∪ B ◦ ))◦ .
Then (67) follows by contraposition.
Next, we show the other inclusion, which is
A ∩ B ⊂ (conv(A◦ ∪ B ◦ ))◦ .

(70)

Suppose that x ∈ A ∩ B = A◦◦ ∩ B ◦◦ . Then for all wA ∈ A◦ , wB ∈ B ◦ , and t ∈ [0, 1], it follows
that
|hx, twA + (1 − t)wB i| ≤ t|hw, wA i| + (1 − t)|hw, wB i| ≤ 1.
Therefore
x ∈ conv(A◦ ∪ B ◦ )◦ .
We have shown that
A ∩ B = conv(A◦ ∪ B ◦ )◦ .
The assertion follows from the definition of polar sets.
By the polar duality in Lemma E.3, the unit ball with respect to the dual of ||| · ||| is given as

 [
 [

conv
conv{ed ⊗ ζn |d ∈ [D]}
Sζ◦n = conv
n∈[N ]

n∈[N ]

= conv{ed ⊗ ζn |n ∈ [N ], d ∈ [D]}.
Collecting the above estimates, we bound the log entropy number in (65) as follows:
log N (B2M D , k · k2∞ , t) . log N (B2M D , ||| · ||| , t/2)


[
. log N conv
Sζ◦n , k · k2 , t/2
n∈[N ]




. log N conv{ed ⊗ ζn |n ∈ [N ], d ∈ [D]}, k · k2 , t/2 ,
where the first inequality follows from (66) and the second inequality holds by the entropy duality
by Artstein et al. [2].
D → `M D so that the standard basis vectors in `N D
Next we define a linear map Q : `N
1
2
1
are mapped to distinct elements in {ed ⊗ ζn | n ∈ [N ], d ∈ [D]}. In this construction, we
only care the one-to-one correspondence and we do not care how we enumerate the elements of

46

KIRYUNG LEE, FELIX KRAHMER, AND JUSTIN ROMBERG

{ed ⊗ ζn | n ∈ [N ], d ∈ [D]}. Although Q is not uniquely determined and there is ambiguity up
D
N D → `M D k = 1.
to a permutation in `N
1 , every map Q constructed as above satisfies that kQ : `1
2
Fix Q and we get
Z
0

1q

log N (B2M D , k

Z

1

r



log N Q(B1N D ), k · k2 , t/2 dt
p0
√ p
. log(N D) log(M D) . D log D log(M D),

· k2,∞ , t)dt .

where the second inequality follows from Corollary A.5.
REFERENCES
[1] A. Ahmed, B. Recht, and J. Romberg, Blind deconvolution using convex programming, IEEE Trans. Inf.
Theory, 60 (2014), pp. 1711–1732.
[2] S. Artstein, V. Milman, and S. J. Szarek, Duality of metric entropy, Annals of mathematics, (2004),
pp. 1313–1328.
[3] L. Balzano and R. Nowak, Blind calibration of sensor networks, in Proc. 6th Int. Conf. Inform. Process. in
Sensor Networks, ACM, 2007, pp. 79–88.
[4] S.-H. Byun, C. M. Verlinden, and K. G. Sabra, Blind deconvolution of shipping sources in an ocean
waveguide, The Journal of the Acoustical Society of America, 141 (2017), pp. 797–807.
[5] B. Carl, Inequalities of Bernstein-Jackson-type and the degree of compactness of operators in banach spaces,
in Annales de l’institut Fourier, vol. 35, 1985, pp. 79–118.
[6] S. Choudhary and U. Mitra, Fundamental limits of blind deconvolution Part I: Ambiguity kernel.
arxiv:1411.3810, November 2014.
[7] S. Choudhary and U. Mitra, Sparse blind deconvolution: What cannot be done, in International Symposium
on Information Theory, IEEE, 2014, pp. 3002–3006.
[8] A. Curtis, P. Gerstoft, H. Sato, R. Snieder, and K. Wapenaar, Seismic interferometryturning noise
into signal, The Leading Edge, 25 (2006), pp. 1082–1092.
[9] K. R. Davidson and S. J. Szarek, Local operator theory, random matrices and Banach spaces, Handbook of
the geometry of Banach spaces, 1 (2001), pp. 317–366.
[10] C. Davis and W. M. Kahan, The rotation of eigenvectors by a perturbation. III, SIAM Journal on Numerical
Analysis, 7 (1970), pp. 1–46.
[11] D. G. Feingold and R. S. Varga, Block diagonally dominant matrices and generalizations of the Gerschgorin
circle theorem, Pacific J. Math, 12 (1962), pp. 1241–1250.
[12] S. Foucart and H. Rauhut, A mathematical introduction to compressive sensing, vol. 1, Birkhäuser Basel,
2013.
[13] J. Garnier and G. Papanicolaou, Passive synthetic aperture imaging, SIAM Journal on Imaging Sciences,
8 (2015), pp. 2683–2705.
[14] G. B. Giannakis and R. W. Heath, Blind indentification of multichannel FIR blurs and perfect image restoration, IEEE Trans. Image Proc., 9 (2000), pp. 1877–1896.
[15] G. H. Golub and C. F. Van Loan, Matrix computations, JHU Press, 2012.
[16] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and
A. Haase, Generalized autocalibrating partially parallel acquisitions (GRAPPA), Magnetic resonance in
medicine, 47 (2002), pp. 1202–1210.
[17] J. H. Gunther and A. L. Swindlehurst, On the use of kernel structure for blind equalization, IEEE Trans.
Sig. Proc., 48 (2000), pp. 799–809.
[18] M. I. Gurelli and C. L. Nikias, EVAM: An eigenvector-based algorithm for multichannel blind deconvolution
of input colored signals, IEEE Trans. Signal Process., 43 (1995), pp. 134–149.
[19] G. Harikumar and Y. Bresler, Perfect blind restoration of images blurred by multiple filters: theory and
efficient algorithms, IEEE Trans. Image Proc., 8 (1999), pp. 202–219.
[20] M. Junge and K. Lee, Generalized notions of sparsity and restricted isometry property. Part I: A unified
framework, arXiv preprint arXiv:1706.09410, (2017).
[21] M. Junge and Q. Zeng, Noncommutative Bennett and Rosenthal inequalities, The Annals of Probability, 41
(2013), pp. 4287–4316.
[22] K. F. Kaaresen and T. Taxt, Multichannel blind deconvolution of seismic signals, Geophysics, 63 (1998),
pp. 2093–2107.

SPECTRAL METHODS FOR PASSIVE IMAGING

47

[23] M. Kech and F. Krahmer, Optimal injectivity conditions for bilinear inverse problems with applications
to identifiability of deconvolution problems, SIAM Journal on Applied Algebra and Geometry, 1 (2017),
pp. 20–37.
[24] V. Koltchinskii, K. Lounici, and A. B. Tsybakov, Nuclear-norm penalization and optimal rates for noisy
low-rank matrix completion, The Annals of Statistics, (2011), pp. 2302–2329.
[25] F. Krahmer, S. Mendelson, and H. Rauhut, Suprema of chaos processes and the restricted isometry property, Comm. Pure Appl. Math., 67 (2014), pp. 1877–1904.
[26] B. Laurent and P. Massart, Adaptive estimation of a quadratic functional by model selection, Annals of
Statistics, (2000), pp. 1302–1338.
[27] K. Lee and M. Junge, RIP-like properties in subsampled blind deconvolution, arXiv preprint arXiv:1511.06146,
(2015).
[28] K. Lee, F. Krahmer, and J. Romberg, An eigen approach to stable multichannel blind deconvolution under
an FIR subspace model, in Sampling Theory and Applications (SampTA), 2017 International Conference
on, 2017.
[29] K. Lee, Y. Li, M. Junge, and Y. Bresler, Blind recovery of sparse signals from subsampled convolution,
IEEE Trans. Inf. Theory, 63 (2017), pp. 802–821.
[30] X. Li, S. Ling, T. Strohmer, and K. Wei, Rapid, robust, and reliable blind deconvolution via nonconvex
optimization, arXiv preprint arXiv:1606.04933, (2016).
[31] Y. Li, K. Lee, and Y. Bresler, Identifiability in blind deconvolution with subspace or sparsity constraints,
IEEE Trans. Inf. Theory, 62 (2016), pp. 4266–4275.
[32] Y. Li, K. Lee, and Y. Bresler, Optimal sample complexity for blind gain and phase calibration, IEEE Trans.
Signal Process., 64 (2016), pp. 5549–5556.
[33] Y. Li, K. Lee, and Y. Bresler, Identifiability and stability in blind deconvolution under minimal assumptions,
IEEE Trans. Inf. Theory, 63 (2017), pp. 4619–4633.
[34] S. Ling and T. Strohmer, Self-calibration via linear least squares, arXiv preprint arXiv:1611.04196, (2016).
[35] H. Liu, G. Xu, L. Tong, and T. Kailath, Recent developments in blind channel equalization : From cyclostationnarity to subspaces, Signal Process., 50 (1996), pp. 83–99.
[36] R. L. Morrison, M. N. Do, and D. C. Munson, Mca: A multichannel approach to sar autofocus, IEEE
transactions on image processing, 18 (2009), pp. 840–853.
[37] E. Moulines, P. Duhamel, J.-F. Cardoso, and S. Mayrargue, Subspace methods for the blind identification
of multichannel FIR filters, IEEE Trans. Signal Process., 43 (1995), pp. 516–525.
[38] H. Q. Nguyen, S. Liu, and M. N. Do, Subspace methods for computational relighting, in Computational
Imaging, 2013, p. 865703.
[39] S. O’Rourke, V. Vu, and K. Wang, Random perturbation of low rank matrices: Improving classical bounds,
arXiv preprint arXiv:1311.2657, (2013).
[40] M. Rudelson and R. Vershynin, Hanson-Wright inequality and sub-gaussian concentration, Electron. Commun. Probab, 18 (2013), pp. 1–9.
[41] K. G. Sabra and D. R. Dowling, Blind deconvolution in ocean waveguides using artificial time reversal, The
Journal of the Acoustical Society of America, 116 (2004), pp. 262–271.
[42] K. G. Sabra, H.-C. Song, and D. R. Dowling, Ray-based blind deconvolution in ocean sound channels, The
Journal of the Acoustical Society of America, 127 (2010), pp. EL42–EL47.
[43] N. Tian, S.-H. Byun, K. Sabra, and J. Romberg, Multichannel myopic deconvolution in underwater acoustic
channels via low-rank recovery, The Journal of the Acoustical Society of America, 141 (2017), pp. 3337–3348.
[44] L. Tong and S. Perreau, Multichannel blind identification: From subspace to maximum likelihood methods,
Proc. IEEE, 86 (1998), pp. 1951–1968.
[45] J. A. Tropp, User-friendly tail bounds for sums of random matrices, Found. Comput. Math., 12 (2012),
pp. 389–434.
[46] V. Vu, Singular vectors under random perturbation, Random Structures & Algorithms, 39 (2011), pp. 526–538.
[47] G. Xu, H. Liu, L. Tong, and T. Kailath, A least-squares approach to blind channel identification, IEEE
Trans. Signal Process., 43 (1995), pp. 2982–2993.
[48] J. Zhou and M. Do, Multidimensional multichannel FIR deconvolution using Gröbner bases, IEEE Trans.
Image Proc., 15 (2006), pp. 2998–3007.
[49] X. Zhu, F. Šroubek, and P. Milanfar, Deconvolving PSFs for a better motion deblurring using multiple
images, Computer Vision–ECCV 2012, (2012), pp. 636–647.

