Solution of linear ill-posed problems using random
dictionaries

arXiv:1605.07913v3 [math.ST] 20 Jun 2017

Pawan Gupta and Marianna Pensky
Department of Mathematics, University of Central Florida
Abstract
In the present paper we consider application of overcomplete dictionaries to solution of general
ill-posed linear inverse problems. In the context of regression problems, there has been enormous
amount of effort to recover an unknown function using such dictionaries. One of the most popular
methods, lasso and its versions, is based on minimizing empirical likelihood and unfortunately,
requires stringent assumptions on the dictionary, the, so called, compatibility conditions. Though
compatibility conditions are hard to satisfy, it is well known that this can be accomplished by using
random dictionaries. In the present paper, we show how one can apply random dictionaries to
solution of ill-posed linear inverse problems. We put a theoretical foundation under the suggested
methodology and study its performance via simulations.
Keywords and phrases: Linear inverse problem; lasso; random dictionaries
AMS (2000) Subject Classification: Primary: 62G05. Secondary: 62C10

1

Introduction

In this paper, we consider solution of a general ill-posed linear inverse problem Qf = q where Q is a
bounded linear operator that does not have a bounded inverse and the right-hand side q is measured
with error. In particular, we consider equation
y = q + ση,

q = Qf ,

(1.1)

where y, q, f , η ∈ Rn , Q ∈ Rn×n . Here, y is observed, q is unobserved, f is the function to be
estimated, σ is the noise level and η is the noise vector which we assume to have the standard
normal distribution. Matrix Q is invertible but its lowest eigenvalue is very small, especially, when n
is relatively large, which makes the problem ill-posed. A general linear inverse problem can usually
be reduced to formulation (1.1) by either expanding y and f over some collection of basis functions
or by measuring them at some set of points.
Solutions of statistical inverse problem (1.1) usually rely on reduction of the problem to the
sequence model by carrying out the singular value decomposition (SVD) (see, e.g., [6], [7], and [15]
and references therein), or its relaxed version, the wavelet-vaguelette decomposition proposed by
Donoho [11] and further studies by Abramovich and Silverman [2]. Another general approach is
Galerkin method with subsequent model selection (see, e.g., [8]).
The advantage of the methodologies listed above is that they are asymptotically optimal in
a minimax sense. The function of interest is usually represented via an orthonormal basis which is
motivated by the form of matrix Q. However, in spite of being minimax optimal in many contexts,
these approaches have drawbacks. In particular, in practical applications, the number of observations
n may be low while noise level σ high. In this situation, if the unknown vector f does not have a
1

relatively compact and accurate representation in the chosen basis, the precision of the resulting
estimator will be poor.
In the last decade, a great deal of effort was spent on recovery of an unknown vector f in
regression setting from its noisy observations using overcomplete dictionaries. In particular, if f has
a sparse representation in some dictionary (a collection of vectors used for the representation of f ),
then f can be recovered with a much better precision than, for example, when it is expanded over
an orthonormal basis. The methodology is based on the idea that the error of an estimator of f
is approximately proportional to the number of dictionary functions that are used for representing
f , therefore, expanding a function of interest over fewer dictionary elements reduces the estimation
error. Similar advantages hold in the case of linear inverse problems (see [14]). However, in order to
represent a variety of functions using a small number of dictionary elements, one needs to consider a
dictionary of much larger size than the number of available observations, the, so called, overcomplete
dictionary.
A variety of techniques have been developed for solution of regression problems using overcomplete dictionaries such as likelihood penalization methods and greedy algorithms. The most popular
of those methods (due to its computational convenience), lasso and its versions, have been used for
solution of a number of theoretical and applied statistical problems (see, e.g., [3], and also [4] and
references therein). However, application of lasso is based on maximizing the likelihood and, unfortunately, relies on stringent assumptions on the dictionary {ϕk }pk=1 , the, so called, compatibility
conditions, for a proof of its optimality. In regression set up, as long as compatibility conditions
hold, lasso identifies a linear combination of the dictionaryp
elements which represent the function of
interest best of all at a ”price” which is proportional to σ n−1 log p where p is the dictionary size
(see, e.g., [4]). Regrettably, while compatibility conditions may be satisfied for the vectors ϕj in the
original dictionary, they usually do not hold for their images Qϕj due to contraction imposed by the
operator Q. Pensky [14] showed how lasso solution can be modified, so that it delivers an optimal
solution, however, compatibility assumptions in [14] remain very complex and hard to verify.
In the recent years it has been discovered that in regression setting, one can satisfy compatibility
conditions for lasso by simply using random dictionaries. In particular, Vershynin [17] provided a
variety of way for construction of such dictionaries, i.e, dictionaries comprised of random vectors. The
methodology of [17], however, is intended for the recovery of a function which is directly observed.
The purpose of the present paper is to explain how random dictionaries can be adopted for solution
of ill-posed linear inverse problems.
To the best of our knowledge, application of random dictionaries to signal recovery has not
been attempted so far since random vectors usually contain “pure noise” and therefore are perceived
as unsuitable for representing a meaningful signal. This is indeed true when one needs to estimate a
simple smooth function which is best represented by a small set of smooth basis functions. However,
when a signal has a more complicated structure, it cannot be expanded over a small number of
basis functions. In this case, a large rich dictionary may be helpful since there is a high chance
that the signal of interest can be well approximated by a linear combination of a small number of
vectors of the dictionary. The advantage of the random dictionaries is that, unlike in the case of
fixed dictionaries, one can work with a dictionary of vary large size which provides a competitive
advantage over standard orthogonal basis based techniques. This benefits can be more significant
when one needs to solve an inverse ill-posed problem since, as it was notes in [14], finding a “low-cost”
representation of a function of interest can significantly improve the accuracy of the solution.
The rest of the paper is organized as follows. Section 2 introduces some notations, formulates
optimization problem with lasso penalty and lists compatibility conditions of [14]. Section 3 contains
the main results of the paper: it explains how one can obtain fast lasso convergence rates by using

2

random dictionaries. Section 4 contains a simulation study which proves that our technique is
competitive. Section 5 concludes the paper with the discussion. Finally, Section 6 contain proofs of
the statements in the paper.

2

Construction of the lasso estimator and a general compatibility
condition

In the paper, we use the following notations.
For any vector t ∈ Rp , denote its ℓ2 , ℓ1 , ℓ0 and ℓ∞ norms by, respectively, ktk, ktk1 , ktk0 and
ktk∞ . For any matix A, denote its ith row and j th column by, Ai· and A·j respectively. Denote
its spectral and Frobenius norms by, respectively, kAk and kAk2 . Denote P = {1, · · · , p}. For any
subset of indices J ⊆ P, subset J c is its complement in P and |J| is its cardinality, so that |P| = p.
Let LJ = Span ϕj , j ∈ J . If J ⊂ P and t ∈ Rp , then tJ ∈ R|J| denotes reduction of vector t to
subset of indices J. Also, ΦJ denotes the reduction of matrix Φ to columns Φ·j with j ∈ J.
Denote by λmin (m; Φ) and λmax (m; Φ) the minimum and the maximum restricted eigenvalues
of matrix ΦT Φ given by
λmin (m; Φ) = minp

t∈R
ktk0 ≤m

tT ΦT Φt
,
ktk22

λmax (m; Φ) = max
p

t∈R
ktk0 ≤m

tT ΦT Φt
.
ktk22

(2.1)

Denote by Φ the dictionary matrix with columns ϕj ∈ Rn , j = 1, · · · , p, where p is possibly
much larger than n and
p
X
tj ϕj = Φt.
(2.2)
ft =
j=1

Let θ be the true vector of coefficients of expansion of f over the dictionary Φ, so that f = Φθ. Let
vectors ψ j be such that QT ψ j = ϕj , where QT is the transpose of matrix Q, and Ψ be a matrix
with columns ψ j , j = 1, · · · , p. Then,
QT Ψ = Φ and

Ψ = Q(QT Q)−1 Φ.

(2.3)

Note that, although f is unknown,
kf − ft k2 = kf k22 + tT ΦT Φt − 2tT ΦT f = kf k22 + tT ΦT Φt − 2tT ΨT Qf

(2.4)

is the sum of the three components where the first one is independent of t, the second one is completely
known, while the last term is of the form 2tT ΨT Qf = 2tT ΨT q and, hence, can be estimated by
2tT ΨT y. Let z be such that
ΨT y = ΦT z.
Therefore, expression tT ΦT Φt−2tT ΨT y is minimized by the same vector t that minimizes kΦt−zk22
where
z = (ΦΦT )−1 ΦΨT y.
(2.5)
Denote νj = kψ j k2 , j = 1, · · · , p, and observe that νj is proportional to the standard deviation
of the j-th component of the vector ΨT y. The value νj can be viewed as a “cost” of using a dictionary
element ϕj in representation of f . Consider a matrix
Υ = diag(ν1 , · · · , νp ) = diag(kψ 1 k2 , · · · , kψ p k2 ).
3

(2.6)

Following [14], we estimate the true vector of coefficients θ as a solution of the quadratic optimization
problem with the weighted lasso penalty

b = arg min kΦt − zk2 + αkΥtk1 .
(2.7)
θ
2
t

Here z is given by (2.5) and α ≥ α0 where

α0 = σ

p

2 n−1 (τ + 1) log p.

(2.8)

Parameter τ > 0 is related to the required probability bound (see formula (6.3) in Section 6 for
b
details). Subsequently, we estimate the unknown solution f by b
f = Φθ.
Note that since we are interested in fθ rather than in the vector θ of coefficients themselves,
we are using lasso for solution of the so called prediction problem where it requires milder conditions
on the dictionary. In fact, it is known (see [14]) that with no additional assumptions, for α ≥ α0 ,
with probability at least 1 − 2p−τ , one has


p
X
νj |tj | .
n−1 kf b − f k22 ≤ inf n−1 kft − f k22 + 4α
(2.9)
θ
t
j=1

It is easy to see that if t = θ, then ft = f . Then, with high probability, the error of the estimator f b
θ
p
P
is proportional to σ n−1 (τ + 1) log p j νj . This is the, so called, slow lasso rate. In order to attain
P
the fast lasso rate proportional to σ 2 n−1 j νj2 , one needs some kind of a compatibility assumption.
Pensky [14] formulated the following compatibility condition: matrices Φ and Υ are such that
for some µ > 1 and any J ⊂ P


dT ΦT Φd · Tr(Υ2J )
2
κ (µ, J) = min d ∈ J (µ, J), kdk2 6= 0 :
> 0.
(2.10)
k(Υd)J k21

where J (µ, J) = d ∈ Rp : k(Υd)J∗c k1 ≤ µk(Υd)J k1 . Pensky [14] proved that, under assumption
(2.10), for α = ̟α0 where ̟ ≥ (µ + 1)/(µ − 1) and α0 is defined in (2.8), with probability at least
1 − 2p−τ , one has


2 K (1 + ̟)2 (τ + 1) log p X
σ
0
νj2  ,
(2.11)
kf b − f k22 ≤ inf kf − fLJ k22 +
θ
J⊆P
κ2 (µ, J)
n
j∈J

where fLJ = projLJ f .
Note, however, that unless matrix Φ has orthonormal columns, assumption (2.10) is hard not
only to satisfy but even to verify since it requires checking it for every subset J in P. Indeed,
sufficient conditions listed in Appendix A1 of [14] rely on the results of Bickel et al. [3] and require
very stringent conditions on λmin (m; Φ) and entries Υ in (2.6). In the present paper, we offer an
alternative to this approach.

3

Lasso solution to linear inverse problems using random dictionaries

An advantage of using random dictionary lies in the fact that one can ensure, with a high probability,
that the dictionary satisfies a restricted isometry condition (see, e.g., [5] or [12]). In particular, if
matrix Φ ∈ Rn×p satisfies the restricted isometry property of order s ≥ 1, then λmin (s; Φ) > 0. The
latter allows one to formulate the following results.
4

Theorem 1 Let θ be the solution of optimization problem (2.7) with α ≥ α0 where α0 is defined in
(2.8). Let Φ ∈ Rn×p be a random dictionary independent of y in (1.1). Denote



X 
(3.1)
J∗ = arg min J ⊂ P : n−1 kf − fLJ k22 + K0 α2
νj2 ,


j∈J

where fLJ = projLJ f and assume that Φ is such that for some s, 1 ≤ s ≤ n/2 and δ, ǫ1 , ǫ2 , ǫ3 ∈ (0, 1),
the following conditions hold
P (λmin (2s; Φ) ≥ 1 − δ) ≥ 1 − ǫ1 ,

P (|J∗ | ≤ s) ≥ 1 − ǫ2 ,


b 0≤s
P kθk
≥ 1 − ǫ3 ,

If K0 ≥ 4/(1 − δ)2 in (3.1), then



X
1
1
P  kf b − f k22 ≤ inf  kf − fLJ k22 + K0 α2
νj2  ≥ 1 − 2p−τ − ǫ1 − ǫ2 − ǫ3 .
J⊆P
n θ
n

(3.2)
(3.3)
(3.4)

(3.5)

j∈J

Note that for α = α0 and K0 = 4/(1 − δ)2 , under assumptions (3.2) – (3.4), formula (3.5) yields
the following result





2
X
1
4σ
1
2 
≥ 1 − 2p−τ − ǫ1 − ǫ2 − ǫ3 . (3.6)
P  kf b − f k22 ≤ inf
kf − fLJ k22 +
ν
j

J⊆P  n
n θ
n (1 − δ)2
j∈J

As Lemma 1 below shows, assumption (3.2) can be guaranteed by choosing a dictionary of a particular
type.

Lemma 1 Let matrix Φ ∈ Rn×p be independent of y and satisfy one of the following conditions:
a) Matrix Φ has independent sub-gaussian isotropic random rows;
b) Matrix Φ has independent sub-gaussian isotropic random columns with unit norms;
√
c) Matrix Φ is obtained as Φ = (c n)−1 DW where W ∈ Rm×p is a matrix with i.i.d. standard
Gaussian entries and columns of the matrix D ∈ Rn×m form a non-random c-tight frame, so that for
any vector x, one has xT DDT x = c2 kxk2 .
If, for some δ ∈ (0, 1) and 1 ≤ s ≤ n/2, one has
n ≥ C1 δ−2 s log(ep/s),

(3.7)

then condition (3.2) holds with ǫ1 ≤ 2 exp(−C2 δ2 n). Here, C1 and C2 depend on the kind of subgaussian variables that are involved in formation of Φ and are independent of n, m, p, s and δ.
Finally, conditions (3.3) and (3.4) can be ensured by restricting the set of solutions t to vectors
with cardinality at most s. In this case, ǫ2 = ǫ3 = 0 and the following corollary of Theorem 1 is valid.
Corollary 1 Let θ be the solution of optimization problem

b = arg min kΦt − zk2 + αkΥtk1 ,
θ
2
t:ktk0 ≤s

5

(3.8)

with α ≥ α0 where α0 is defined in (2.8). Let Φ ∈ Rn×p be one of the random random dictionaries
defined in Lemma 1. If, for some δ ∈ (0, 1), condition (3.7) holds, then



2
X
1
1
4α
P  kf b − f k22 ≤ inf  kf − fLJ k22 +
νj2  ≥ 1 − 2p−τ − 2 exp(−C2 δ2 n),
(3.9)
J ⊆P
n θ
n
(1 − δ)2
j∈J

|J|≤s

where C2 depends on the kind of sub-gaussian variables that are involved in formation of Φ and is
independent of n, m, p, s and δ.
Note that case c) above offers a structured random dictionary since each of its elements is a
linear combination of smooth functions.

4

Simulation studies

In order to evaluate the performance of the procedure suggested in this paper we carried out a limited
simulation study. For our study, we chose three sample sizes n = 32, n = 64 and n = 128. We first
generated a true vector f using MakeSignal program in the package Wavelab 850. We then generated
the matrix Q in (1.1) as Q = UΛUT where
√ U is an (n × n) random orthogonal matrix and Λ is
a diagonal matrix with entries Λii = 1/ i, i = 1, 2, · · · , n. Using Q we obtained the unobserved
vector q as
q = Qf .
At last, for generating the data y we added Gaussian random noise to q. For this purpose, we chose
particular values of the Signal to Noise Ratio (SNR) and obtained σ as the ratio of the standard
deviation of q and the SNR. Vector y was then calculated at n observation points as y = q + σ η
where η ∈ Rn is a standard normal vector. Finally, we ran simulations for two noise levels: SNR =
3 and SNR = 5.
We compared the estimators of f based on random dictionaries with the estimator of f based on
the Singular Value Decomposition (SVD). For our simulations we have created three different n × p
random dictionaries with p = 5000: (a) two purely random dictionaries with, respectively, the i.i.d.
standard Gaussian entries and the i.i.d. sparse Bernoulli entries; (b) the fusion of the fixed dictionary
and the random dictionary that follows case c) in Lemma 1 with D being the Haar dictionary. The
sparse Bernoulli variable is defined as
 q
3
1


 − n with probability 6
0q
with probability 23
X=


3

with probability 1
n

6

For creating the fusion dictionary, we first generated the orthogonal matrix of the Haar wavelet
transform D using MakeWavelet function, so that m = n and c = 1. Then we obtained the dictionary
Φ following part c) of the Lemma 1 using the n × p matrix W with the i.i.d. normal entries.
We obtained matrix Ψ of the inverse images as the numerical solution of the exact equation
QT Ψ = Φ and calculated vector z with elements (2.5). For the sake of obtaining a solution of
optimization problem (2.7), we used function LassoWeighted in SPAMS MatLab toolbox (see [13]).
In order to evaluate the value of the lasso parameter α, we calculated αmax as the value of
α that guarantees that all coefficients in the model vanish. We created a grid of the values of
αk = αmax ∗ k/N , k = 1, · · · , N , with N = 200. As a result, we obtained a collection of estimators
6

6

3

3

2.5

5
2

2

4
3

1.5

1

1
2
0.5

0
1

0
0

-1

-0.5

-1

-1

-2
-2

-1.5

-3

-3
0

10

20

30

40

50

60

70

-2
0

10

20

30

40

50

60

70

0

10

20

30

40

50

60

70

Figure 1: Test signals WernerSorrows (left), MishMash (middle) and Chirps (right) with n = 64.
b = θ(α
b k ). For the purpose of choosing the most appropriate value of k, we estimated α as α̂ = αb in
θ
k
two ways: one using the oracle value of α and another using the estimated value of α. We found oracle
koracle /N using the value b
koracle that guarantees the most accurate
value of α as αoracle = αmax ∗ b
estimator of f :
b
b k )k2 .
koracle = arg min kf − Φθ(α
k

Since the vector f is unavailable in real life, we find the estimated value α̂est = αmax ∗ b
kest /N of α
using


1
2
2 −1
b
ky − q̂(αk )k2 + 2σ n p̂k ,
kest = arg min
k
n

b k ) is the estimator of q based on the lasso estimator obtained with the
where q̂(αk ) = QΦθ(α
b k ).
parameter αk and p̂k is the number of nonzero components of θ(α
We compared the estimators f̂RN , f̂RB , f̂RH of f based, respectively, on Gaussian, Bernoulli
and Haar fusion random dictionaries described above with f̂SV D , the estimator based on the singular
value decomposition (SVD). Initially we considered wavelet estimator of f using Daubechies wavelet
of order 8, but we discarded it due to its poor performance with respect to the estimators considered
for comparison. For finding f̂SV D , we used the oracle number Koracle of eigenbasis functions. We
obtained Koracle as the number of eigenbasis functions that minimizes the difference between f̂SV D
and the true function f which is unavailable in a real life setting.
Table 1 below compares the accuracies of the estimators based on random dictionaries with the
SVD estimator. Precision of an estimator f̂ is measured by R(f̂ ) = n−1/2 kf̂ − f k2 , the estimated L2 norm of the difference between the estimator f̂ and the true vector f averaged over 50 simulation runs
(with the standard deviations listed in parentheses). For all the three estimator based on random
oracle , f̂ oracle ,
dictionaries, we report the errors with both the oracle and the estimated values of α, f̂RN
RB
oracle and f̂ est , f̂ est , f̂ est , respectively. We carried out simulations with three types of test functions
f̂RH
RN RB RH
WernerSorrows, MishMash and Chirps. The test signals are presented in Figure 1.
From Table 1 it follows that all the random dictionary based estimators are more accurate than
oracle , f̂ oracle and f̂ oracle over f̂
the SVD estimator. The advantage of f̂RN
SV D is more significant than
RB
RH
est
est
est
that of f̂RN , f̂RB and f̂RH since the latter estimators loose accuracy because of suboptimal choices
of the parameter α. Nevertheless, in majority of cases, they still exhibit better precision than f̂SV D
although this is not entirely fair comparison since f̂SV D is based on the oracle choice of parameter
K. This is due to the fact that large random dictionaries provide a more sparse representation of f .

7

0.5127
(0.0767)
0.3956
(0.0632)
0.6074
( 0.0541)
0.3208
(0.0350)
0.6072
(0.0431 )
0.3807
(0.0267)

W ernerSorrows
oracle
est
f̂RB
f̂RB
0.4910
0.5125
(0.0736) (0.0803)
0.3810
0.3966
(0.0654) (0.0637)
0.5874
0.6112
(0.0474) (0.0537)
0.3166
0.3210
(0.0360) (0.0364)
0.5780
0.6082
(0.0394) (0.0442)
0.3717
0.3801
(0.0244) (0.0268)

oracle
f̂RH
0.4837
(0.0724)
0.3801
(0.0622)
0.5759
(0.0461)
0.3160
(0.0349)
0.5744
(0.0393)
0.3717
(0.0250)

est
f̂RH
0.4981
(0.0716)
0.3917
(0.0631)
0.5929
(0.0510)
0.3199
( 0.0359)
0.5901
(0.0417)
0.3757
(0.0259)

f̂SV D
0.5155
(0.0721)
0.4020
(0.0632)
0.6114
(0.0577)
0.3254
(0.0359)
0.6218
(0.0440)
0.3850
(0.0271)

0.5409
(0.0811)
0.3523
(0.0588)
0.6131
(0.0520)
0.3039
(0.0298)
0.5112
(0.0393)
0.3385
(0.0244)

0.5628
(0.0871)
0.3595
(0.0578)
0.6400
(0.0617)
0.3086
(0.0315)
0.5252
(0.0409)
0.3430
(0.0240)

M ishM ash
0.5432
0.5631
(0.0800) (0.0854)
0.3537
0.3613
(0.0561) (0.0566)
0.6158
0.6391
(0.0552) (0.0595)
0.3039
0.3083
(0.0283) (0.0300)
0.5106
0.5246
(0.0381) (0.0412)
0.3380
0.3429
(0.0241) (0.0242)

0.5406
(0.0848)
0.3532
( 0.0554)
0.6145
(0.0541)
0.3024
(0.0293)
0.5110
(0.0394)
0.3383
(0.0237)

0.5562
(0.0846)
0.3592
(0.0560)
0.6352
(0.0606)
0.3067
(0.0294)
0.5185
(0.0402)
0.3410
( 0.0235)

0.5977
(0.0874)
0.3657
(0.0549)
0.6599
(0.0621)
0.3107
(0.0297)
0.5410
(0.0415)
0.3460
(0.0245)

0.4300
(0.0630)
0.2872
(0.0461)
0.3979
(0.0369)
0.2735
(0.0314)
0.3878
(0.0288)
0.2428
(0.0182)

0.4430
(0.0646)
0.2977
(0.0459)
0.4117
(0.0404)
0.2782
(0.0345)
0.3955
( 0.0310)
0.2460
(0.0179)

Chirps
0.4320
0.4431
(0.0602) (0.0598)
0.2880
0.2960
(0.0448) ( 0.0448)
0.3983
0.4121
(0.0393) (0.0414)
0.2733
0.2778
(0.0327) (0.0351)
0.3867
0.3953
(0.0282) (0.0289)
0.2431
0.2456
(0.0180) (0.0179)

0.4307
(0.0613)
0.2871
(0.0457)
0.3986
(0.0406)
0.2723
(0.0320)
0.3868
( 0.0289)
0.2428
(0.0180)

0.4443
(0.0623)
0.2951
(0.0444)
0.4066
(0.0404)
0.2767
(0.0341)
0.3953
( 0.0285)
0.2460
(0.0181)

0.4607
(0.0607)
0.3018
(0.0443)
0.4297
(0.0391)
0.2789
( 0.0347)
0.4069
(0.0293)
0.2473
(0.0173)

oracle
f̂RN

est
f̂RN

n = 32,
SN R = 3
n = 32,
SN R = 5
n = 64,
SN R = 3
n = 64,
SN R = 5
n = 128,
SN R = 3
n = 128,
SN R = 5

0.4910
(0.0715)
0.3838
(0.0642)
0.5818
(0.0511)
0.3152
(0.0351)
0.5761
(0.0406)
0.3730
(0.0251)

n = 32,
SN R = 3
n = 32,
SN R = 5
n = 64,
SN R = 3
n = 64,
SN R = 5
n = 128,
SN R = 3
n = 128,
SN R = 5

n = 32,
SN R = 3
n = 32,
SN R = 5
n = 64,
SN R = 3
n = 64,
SN R = 5
n = 128,
SN R = 3
n = 128,
SN R = 5

Table 1: The average values of the errors R(f̂ ) evaluated over 50 simulation runs of the estimators
for various test signals (with the standard deviations of the errors listed in the parentheses).
8

5

Discussion

In the present paper we provided a new approach for the solution of a general ill-posed linear inverse
problem. The underlying idea is to use lasso technique for estimating the function of interest by
representing it as a sparse linear combination of elements of a random overcomplete dictionary. The
advantage of choosing a random dictionary over any other overcomplete dictionary is that one can
construct it in such a way that it satisfies restricted isometry condition with a high probability and,
therefore, ensures that the compatibility condition (which guarantees fast convergence rates for lasso)
also holds.
We provide theoretical justification for application of the lasso technique with the random
dictionaries for solution of the linear inverse problems. We also support our theory by the simulation
studies which show that the proposed estimators have higher accuracy than the SVD estimators in
spite of the fact that the SVD estimators are based on the oracle parameter choices. For this reason,
the advantage of the random dictionary based estimators is more significant when they are likewise
constructed with the oracle choices of parameter α. In fact, this is the part where our method has
some room for improvement: since our procedure for estimating parameter α is rather elementary, it
can be fine-tuned using more advanced techniques.

Acknowledgements
Marianna Pensky and Pawan Gupta were partially supported by National Science Foundation (NSF),
grants DMS-1106564 and DMS-1407475.

6

Proofs

Proof of Theorem 1. The beginning of the proof is similar to the proof of Lemma 2 in [14].
However, for completeness, we provide the complete proof here.
Let θ be the true parameter vector, so that f = fθ = Φθ. Denote ζ = ΨT η. Then, it is easy
to check that
ΦT (z − f ) = ΨT (y − Qf ) = σζ.
Following [10], by K-K-T condition, we derive that for any t ∈ Rp
bT ΦT (z − Φθ)
b = α
θ
b ≤ α
tT ΦT (z − Φθ)

p
X

j=1
p
X
j=1

νj |θbj |
νj |tj |,

so that, subtracting the first line from the second, we obtain
b − Φt)T (Φθ
b − z) ≤ α
(Φθ

p
X
j=1

νj (|tj | − |θbj |).

(6.1)

b −Φt)T (Φθ
b −Φθ) ≤ σ(θ
b −t)T ζ +α Pp νj (|tj |−|θbj |). Since for any u, v ∈ Rp
Then, (6.1) yields (Φθ
j=1


b − Φt and u = Φθ
b − Φθ for any t ∈ Rp
one has vT u = 12 kvk2 + kuk2 − kv − uk2 , choosing v = Φθ
9

obtain
b − t)k2 ≤ kft − f k2 + 2σ(θ
b − t)T ζ + 2α
kf b − f k + kΦ(θ
θ
2

p
X
j=1

νj (|tj | − |θbj |).

(6.2)

By definition of ζ, for any j = 1, · · · , p, one has ζj ∼ N (0, νj2 ). Hence, on the set
Ω0 =



ω:

max (ν −1 |ζj |)
1≤j≤p j

≤

p

2(τ + 1) log p



with

P(Ω0 ) ≥ 1 − 2p−τ

(6.3)

p
b − t)T ζ| ≤ 2(τ + 1) log p Pp νj |θbj − tj | = α0 Pp νj |θbj − tj |. Combining the last
one obtains |(θ
j=1
j=1
inequality with (6.2) obtain that, for any α > 0, on the set Ω0 ,
b − t)k2 ≤ kft − f k2 + 2α
kf b − f k + kΦ(θ
θ
2

p
X
j=1

νj (|tj | − |θbj |) + 2α0

p
X
j=1

νj |θbj − tj |.

(6.4)

n
o
b 0 ≤ s . Choose
Denote Ω1 = {ω : λmin (2s; Φ) ≥ 1 − δ}, Ω2 = {ω : |J∗ | ≤ s} and Ω3 = ω : kθk
t such that ft = projLJ∗ f = fLJ∗ and note that tj = 0 for j ∈ J∗c . Then, due to α ≥ α0 and
||θbj − tj | ≤ |θbj | + |tj |, obtain
b − t)k2 ≤ kft − fL k2 + 4α
kf b − f k2 + kΦ(θ
J∗
θ

X

j∈J∗

νj |θbj − tj |.

(6.5)

Consider the set Ω = Ω0 ∩ Ω1 ∩ Ω2 ∩ Ω3 and note that P(Ω) ≥ 1 − 2p−τ − ǫ1 − ǫ2 − ǫ3 . If
b − tk0 ≤ 2s and, hence,
ω ∈ Ω, then kθ
4α

X

j∈J∗



νj |θbj − tj | ≤ 4α 

X

j∈J∗

1/2

νj2 

2
X
b − t)k
kΦ(θ
b − t)k2 + 4α
νj2 .
≤ kΦ(θ
2
λmin (2s; Φ)
(1 − δ)
j∈J∗

Plugging the last inequality into (6.5) and recalling the definition of J∗ , we derive (3.5).
Proof of Lemma 1. In cases a) and b), λmin (m; Φ) ≥ 1 − δ is ensured by Theorem 5.65 of
Vershynin [17]. In case c), note that entries of matrix Φ are uncorrelated and, hence, are independent
Gaussian variables due to
Cov (Φik · Φjl ) =

m
m
1 X X
Dir1 Djr2 I(r1 = r2 )I(k = l) = I(i = j)I(k = l).
c2
r1 =1 r2 =1

Moreover, matrix Φ has isotropic rows since
m
m
1 X X
Cov (Φih · Φjl ) = 2
Dir1 Djr2 I(r1 = r2 )I(h = l) = I(i = j)I(h = l).
c
r1 =1 r2 =1

Therefore, λmin (m; Φ) ≥ 1 − δ by Theorem 5.65 of [17].

10

References
[1] Abramovich, F., Pensky, M., Rozenholc, Y. (2013) Laplace deconvolution with noisy observations. Electronic Journal of Statistics, 7, 1094-1128
[2] Abramovich, F., Silverman, B. W. (1998). Wavelet decomposition approaches to statistical inverse problems. Biometrika , 85, 115–129.
[3] Bickel, P.J., Ritov, Y., Tsybakov, A. (2009) Simultaneous analysis of Lasso and Dantzig selector.
Ann. Statist., 37, 1705 - 1732.
[4] Bühlmann, P., van de Geer, S. (2011) Statistics for High-Dimensional Data: Methods, Theory
and Applications. Springer.
[5] Candès, E. J., Eldar, Y., Needell, D., Randall, P. (2010) Compressed sensing with coherent and
redundant dictionaries. Appl. Computat. Harmonic Anal., 31, 59–73.
[6] Cavalier, L., Golubev, G.K., Picard, D., Tsybakov, A.B. (2002) Oracle inequalities for inverse
problems. Ann. Statist., 30, 843-874.
[7] Cavalier, L., Reiss, M. (2014) Sparse model selection under heterogeneous noise: Exact penalisation and data-driven thresholding. Electronic Journ. Statist., 8, 432-455.
[8] Cohen, A., Hoffmann, M., Reiss, M. (2004) Adaptive wavelet Galerkin methods for linear inverse
problems. SIAM Journ. Numer. Anal., 42, 1479–1501.
[9] Comte, F., Cuenod, C.-A., Pensky, M., Rozenholc, Y. (2017) Laplace deconvolution on the basis
of time domain data and its application to Dynamic Contrast Enhanced imaging. Journ. Royal
Stat. Soc., Ser.B, 79, 6994.
[10] Dalalyan, A.S., Hebiri, M., Lederer, J. (2014) On the prediction performance of the Lasso.
Bernoulli, 23, 552-581.
[11] Donoho, D.L. (1995). Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition. Applied and Computational Harmonic Analysis, 2, 101–126.
[12] Foucart, S., Rauhut, H. (2013) A Mathematical Introduction to Compressive Sensing. Springer,
New York.
[13] Mairal,
J. (2014) SPAMS: a Sparse
http://spams-devel.gforge.inria.fr

Modeling

Software,

MatLab

toolbox.

[14] Pensky, M. (2016) Solution of linear ill-posed problems using overcomplete dictionaries. Ann.
Statist., 44, 1739-1764.
[15] Tropp, J.A., Wright, S. J. (2010) Computational methods for sparse solution of linear inverse
problems. Proc. IEEE, special issue, ”Applications of sparse representation and compressive
sensing”, 98, 948-958.
[16] Vareschi T. (2013) Noisy Laplace deconvolution with error in the operator. J. Statist. Plan. Inf.,
157-158, 16-35.
[17] Vershynin, R. (2012) Introduction to the non-asymptotic analysis of random matrices. In Compressed Sensing, Theory and Applications, ed. Y. Eldar and G. Kutyniok, Chapter 5. Cambridge
University Press.
11

