Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This tour explores theoritical garantees for the performance of recovery using $\ell^1$ minimization.
using PyPlot
using NtToolBox
# using Autoreload
# arequire("NtToolBox")
WARNING: Method definition ndgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:3 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:3. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:6 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:6. WARNING: Method definition ndgrid_fill(Any, Any, Any, Any) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:13 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:13. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}...) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:19 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:19. WARNING: Method definition meshgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:33 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:33. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:36 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:36. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:44 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:44.
We consider the inverse problem of estimating an unknown signal $x_0 \in \RR^N$ from noisy measurements $y=\Phi x_0 + w \in \RR^P$ where $\Phi \in \RR^{P \times N}$ is a measurement matrix with $P \leq N$, and $w$ is some noise.
This tour is focused on recovery using $\ell^1$ minimization $$ x^\star \in \uargmin{x \in \RR^N} \frac{1}{2}\norm{y-\Phi x}^2 + \la \normu{x}. $$
Where there is no noise, we consider the problem $ \Pp(y) $ $$ x^\star \in \uargmin{\Phi x = y} \normu{x}. $$
We are not concerned here about the actual way to solve this convex problem (see the other numerical tours on sparse regularization) but rather on the theoritical analysis of wether $x^\star$ is close to $x_0$.
More precisely, we consider the following three key properties
\Pp(y) $ for $y=\Phi x_0$.
O(\norm{w})$ for $y=\Phi x_0+w$ if $\norm{w}$ is smaller than an arbitrary small constant that depends on $x_0$ if $\la$ is well chosen according to $\norm{w}$.
arbitrary.
Note that noise robustness implies identifiability, but the converse is not true in general.
The simplest criteria for identifiality are based on the coherence of the matrix $\Phi$ and depends only on the sparsity $\norm{x_0}_0$ of the original signal. This criteria is thus not very precise and gives very pessimistic bounds.
The coherence of the matrix $\Phi = ( \phi_i )_{i=1}^N \in \RR^{P \times N}$ with unit norm colum $\norm{\phi_i}=1$ is $$ \mu(\Phi) = \umax{i \neq j} \abs{\dotp{\phi_i}{\phi_j}}. $$
Compute the correlation matrix (remove the diagonal of 1's).
remove_diag = C -> C - diagm(diag(C))
Correlation = Phi -> remove_diag(abs(Phi'*Phi))
(::#3) (generic function with 1 method)
Compute the coherence $\mu(\Phi)$.
mu = Phi -> maximum(Correlation(Phi))
(::#5) (generic function with 1 method)
The condition $$ \normz{x_0} < \frac{1}{2}\pa{1 + \frac{1}{\mu(\Phi)}} $$ implies that $x_0$ is identifiable, and also implies to robustess to small and bounded noise.
Equivalently, this condition can be written as $\text{Coh}(\normz{x_0})<1$ where $$ \text{Coh}(k) = \frac{k \mu(\Phi)}{ 1 - (k-1)\mu(\Phi) } $$
Coh = (Phi, k) -> (k*mu(Phi))/(1 - (k - 1)*mu(Phi))
(::#7) (generic function with 1 method)
Generate a matrix with random unit columns in $\RR^P$.
normalize = Phi -> Phi ./ repeat(sqrt(sum(Phi.^2, 1)), inner = [size(Phi)[1], 1]);
PhiRand = (P, N) -> normalize(randn(P, N))
Phi = PhiRand(250, 1000)
# P = 250; N = 1000
# Phi = randn(P, N)
# repeat(sqrt(sum(Phi.^2, 1)), inner = [size(Phi)[1], 1])
# sum(Phi.^2, 1)
250×1000 Array{Float64,2}: -0.00153558 -0.0873719 -0.0267141 … -0.00418114 0.0431689 0.016081 0.0658395 0.0183717 -0.0977312 -0.173786 -0.0202818 0.10513 -0.0061723 0.0903024 0.0475056 0.0174428 0.0767283 0.0283078 0.0351142 0.0178382 -0.0156658 -0.113847 0.0829878 -0.0569733 -0.082541 -0.159691 0.0178282 0.102784 … 0.0324298 0.0500197 0.109889 -0.00281318 0.0497717 0.0133811 0.0990111 0.158323 -0.00278128 0.0263971 -0.00234905 0.0866833 0.000221441 0.00898388 0.0364112 -0.0818329 -0.0140391 -0.023077 -0.101275 -0.0164376 -0.0694768 0.056435 0.0391143 0.0801067 -0.00795326 … 0.0105741 -0.0107431 0.0457144 -0.00871715 -0.0343846 0.0428562 -0.0408472 0.00232064 0.0631691 -0.0595921 -0.000329357 0.0098766 ⋮ ⋱ 0.0489705 -0.0196438 0.00785909 1.66286e-5 -0.0416976 0.0460609 -0.0325308 -0.108784 -0.0299601 -0.0500156 0.0730461 -0.0729764 -0.0382131 … 0.00106382 -0.0420094 0.0996469 -0.0980691 -0.0351391 0.0462986 0.0486122 0.0918515 0.106906 -0.0167351 -0.0218581 -0.101535 -0.0554057 0.00990167 0.0265645 -0.0473953 0.0306916 0.0601045 -0.0544066 -0.0442649 -0.0428824 0.0270975 0.0339008 0.0220354 0.0601455 … -0.0811251 -0.013028 0.0160892 0.0403653 0.00898637 0.0820708 0.00224229 -0.0368455 -0.040219 -0.0506181 0.0362769 -0.0364124 0.0991768 0.0148596 0.067233 0.0141599 -0.0624272 0.00481869 -0.0303919 -0.00353946 0.0225215 -0.056659
Compute the coherence and the maximum possible sparsity to ensure recovery using the coherence bound.
c = mu(Phi)
println(@sprintf("Coherence: %.2f", c))
println(@sprintf("Sparsity max: %d", floor(1/2*(1 + 1/c))))
Coherence: 0.32 Sparsity max: 2
Exercise 1
Display how the average coherence of a random matrix decays with the redundancy $\eta = P/N$ of the matrix $\Phi$. Can you derive an empirical law between $P$ and the maximal sparsity?
include("NtSolutions/sparsity_6_l1_recovery/exo1.jl")
## Insert your code here.
In the following we will consider the support $$ \text{supp}(x_0) = \enscond{i}{x_0(i) \neq 0} $$ of the vector $x_0$. The co-support is its complementary $I^c$.
supp = s -> find(abs(s) .> 1e-5)
cosupp = s -> find(abs(s) .< 1e-5);
Given some support $ I \subset \{0,\ldots,N-1\} $, we will denote as $ \Phi = (\phi_m)_{m \in I} \in \RR^{N \times \abs{I}}$ the sub-matrix extracted from $\Phi$ using the columns indexed by $I$.
J.J. Fuchs introduces a criteria $F$ for identifiability that depends on the sign of $x_0$.
J.J. Fuchs. Recovery of exact sparse representations in the presence of bounded noise. IEEE Trans. Inform. Theory, 51(10), p. 3601-3608, 2005
Under the condition that $\Phi_I$ has full rank, the $F$ measure of a sign vector $s \in \{+1,0,-1\}^N$ with $\text{supp}(s)=I$ reads $$ \text{F}(s) = \norm{ \Psi_I s_I }_\infty \qwhereq \Psi_I = \Phi_{I^c}^* \Phi_I^{+,*} $$ where $ A^+ = (A^* A)^{-1} A^* $ is the pseudo inverse of a matrix $A$.
The condition $$ \text{F}(\text{sign}(x_0))<1 $$ implies that $x_0$ is identifiable, and also implies to robustess to small noise. It does not however imply robustess to a bounded noise.
Compute $\Psi_I$ matrix.
PsiI = (Phi,I) -> Phi[:, setdiff(1:size(Phi)[2], I)]' * pinv(Phi[:,I])';
Compute $\text{F}(s)$.
F = (Phi, s) -> norm(PsiI(Phi, supp(s))*s[supp(s)], Inf)
(::#19) (generic function with 1 method)
The Exact Recovery Criterion (ERC) of a support $I$, introduced by Tropp in
J. A. Tropp. Just relax: Convex programming methods for identifying sparse signals. IEEE Trans. Inform. Theory, vol. 52, num. 3, pp. 1030-1051, Mar. 2006.
Under the condition that $\Phi_I$ has full rank, this condition reads $$ \text{ERC}(I) = \norm{\Psi_{I}}_{\infty,\infty} = \umax{j \in I^c} \norm{ \Phi_I^+ \phi_j }_1. $$ where $\norm{A}_{\infty,\infty}$ is the $\ell^\infty-\ell^\infty$ operator norm of a matrix $A$.
erc =(Phi, I) -> norm(PsiI(Phi,I), Inf);
The condition $$ \text{ERC}(\text{supp}(x_0))<1 $$ implies that $x_0$ is identifiable, and also implies to robustess to small and bounded noise.
One can prove that the ERC is the maximum of the F criterion for all signs of the given support $$ \text{ERC}(I) = \umax{ s, \text{supp}(s) \subset I } \text{F}(s). $$
The weak-ERC is an approximation of the ERC using only the correlation matrix $$ \text{w-ERC}(I) = \frac{ \umax{j \in I^c} \sum_{i \in I} \abs{\dotp{\phi_i}{\phi_j}} }{ 1-\umax{j \in I} \sum_{i \neq j \in I} \abs{\dotp{\phi_i}{\phi_j}} }$$
g = (C, I) -> sum(C[:,I], 2)
werc_g = (g, I, J) -> maximum(g[J]) / (1 - maximum(g[I]))
werc = (Phi, I) -> werc_g(g(Correlation(Phi), I), I, setdiff(1:size(Phi)[2], I) )
(::#27) (generic function with 1 method)
One has, if $\text{w-ERC}(I)>0$, for $ I = \text{supp}(s) $, $$ \text{F}(s) \leq \text{ERC}(I) \leq \text{w-ERC}(I) \leq \text{Coh}(\abs{I}). $$
This shows in particular that the condition $$ \text{w-ERC}(\text{supp}(x_0))<1 $$ implies identifiability and robustess to small and bounded noise.
Exercise 2
Show that this inequality holds on a given matrix. What can you conclude about the sharpness of these criteria ?
include("NtSolutions/sparsity_6_l1_recovery/exo2.jl")
N = 2000, P = 1990, |I| = 6 F(s) = 0.19 ERC(I) = 0.23 w-ERC(s) = 0.28 Coh(|s|) = 1.79
## Insert your code here.
N = 2000
P = N - 10
Phi = PhiRand(N, P)
s = zeros(N)
s[1:6] = 1
I = supp(s)
k = length(I)
println(@sprintf("N = %d, P = %d, |I| = %d", N, P, k))
println(@sprintf("F(s) = %.2f", F(Phi, s)))
println(@sprintf("ERC(I) = %.2f", erc(Phi, I)))
println(@sprintf("w-ERC(s) = %.2f", werc(Phi, I)))
println(@sprintf("Coh(|s|) = %.2f", Coh(Phi, k)))
N = 2000, P = 1990, |I| = 6 F(s) = 0.20 ERC(I) = 0.23 w-ERC(s) = 0.26 Coh(|s|) = 1.54
Exercise 3
For a given matrix $\Phi$ generated using PhiRand, draw as a function of the sparsity $k$ the probability that a random sign vector $s$ of sparsity $\norm{s}_0=k$ satisfies the conditions $\text{F}(x_0)<1$, $\text{ERC}(x_0)<1$ and $\text{w-ERC}(x_0)<1$
include("NtSolutions/sparsity_6_l1_recovery/exo3.jl")
## Insert your code here.
N = 600
P = Int(N/2)
Phi = PhiRand(P, N)
klist = Array{Int64,1}(round(linspace(1, P/7., 20)))
ntrials = 60
proba = zeros(length(klist), 3)
for i in 1:length(klist)
proba[i, 1:3] = 0
k = Int(klist[i])
for j in 1:ntrials
s = zeros(N)
I = randperm(N)
I = I[1:k]
l = randn(k, 1)
s[I] = l./abs(l)
proba[i, 1] = proba[i, 1] + (F(Phi, s) .< 1)
proba[i, 2] = proba[i, 2] + (erc(Phi, I) .< 1)
proba[i, 3] = proba[i, 3] + (werc(Phi, I) .> 0).*(werc(Phi, I) .< 1)
end
end
figure(figsize = (8, 5))
plot(klist, proba/ntrials, linewidth = 2)
xlabel("k")
legend(["F <1", "ERC <1", "w-ERC <1"])
title(@sprintf("N = %d, P = %d", N, P))
show()
The restricted isometry constants $\de_k^1,\de_k^2$ of a matrix $\Phi$ are the smallest $\de^1,\de^2$ that satisfy $$ \forall x \in \RR^N, \quad \norm{x}_0 \leq k \qarrq (1-\de^1)\norm{x}^2 \leq \norm{\Phi x}^2 \leq (1+\de^2)\norm{x}^2. $$
E. Candes shows in
E. J. Cand s. The restricted isometry property and its implications for compressed sensing. Compte Rendus de l'Academie des Sciences, Paris, Serie I, 346 589-592
that if $$ \de_{2k} \leq \sqrt{2}-1 ,$$ then $\norm{x_0} \leq k$ implies identifiability as well as robustness to small and bounded noise.
The stability constant $\la^1(A), \la^2(A)$ of a matrix $A = \Phi_I$ extracted from $\Phi$ is the smallest $\tilde \la_1,\tilde \la_2$ such that $$ \forall \al \in \RR^{\abs{I}}, \quad (1-\tilde\la_1)\norm{\al}^2 \leq \norm{A \al}^2 \leq (1+\tilde \la_2)\norm{\al}^2. $$
These constants $\la^1(A), \la^2(A)$ are easily computed from the largest and smallest eigenvalues of $A^* A \in \RR^{\abs{I} \times \abs{I}}$
minmax = v -> (1 - minimum(v), maximum(v) - 1)
ric = A -> minmax(eig(A'*A)[1])
(::#31) (generic function with 1 method)
The restricted isometry constant of $\Phi$ are computed as the largest stability constants of extracted matrices $$ \de^\ell_k = \umax{ \abs{I}=k } \la^\ell( \Phi_I ). $$
The eigenvalues of $\Phi$ are essentially contained in the interval $ [a,b] $ where $a=(1-\sqrt{\be})^2$ and $b=(1+\sqrt{\be})^2$ with $\beta = k/P$ More precisely, as $k=\be P$ tends to infinity, the distribution of the eigenvalues tends to the Marcenko-Pastur law $ f_\be(\la) = \frac{1}{2\pi \be \la}\sqrt{ (\la-b)^+ (a-\la)^+ }. $
Exercise 4
Display, for an increasing value of $k$ the histogram of repartition of the eigenvalues $A^* A$ where $A$ is a Gaussian matrix of size $(P,k)$ and variance $1/P$. For this, accumulate the eigenvalues for many realizations of $A$.
include("NtSolutions/sparsity_6_l1_recovery/exo4.jl")
# klist = [10, 30, 50]
# P = 200
# ntrials = 200
# tmin = 0
# tmax = 2.5
# q = 50
# t = linspace(tmin, tmax, q)
# t1 = linspace(tmin, tmax, 1000)
# dt = (tmax - tmin)/q
# for j in 1 : length(klist)
# k = klist[j]
# # simulation
# v = []
# for i in 1 : ntrials
# v = [v; svd(randn(P, k)./sqrt(P))[2].^2]
# end
# figure(figsize = (10, 10))
# subplot(length(klist), 1, j)
# h = hist(v, t)[2]
# h = h/sum(h)/dt
# h = [h; 0]
# bar(t[1 : end], h, width = 1/20, color = "darkblue", edgecolor = "black")
# # theoritical law
# beta = k/P
# a = (1 - sqrt(beta))^2
# b = (1 + sqrt(beta))^2
# z = sqrt(max(t1 - a, zeros(length(t1))).*max(b - t1, zeros(length(t1))))./(2*pi.*beta.*t1)
# plot(t1, z, "r", linewidth = 3)
# xlim(tmin, tmax)
# ylim(0, maximum(h)*1.05)
# title(@sprintf("P = %d, k = %d", P, k))
# show()
# end
WARNING: hist(...) and hist!(...) are deprecated. Use fit(Histogram,...) in StatsBase.jl instead. in depwarn(::String, ::Symbol) at ./deprecated.jl:64 in #hist!#994(::Bool, ::Function, ::Array{Int64,1}, ::Array{Any,1}, ::LinSpace{Float64}) at ./deprecated.jl:629 in hist(::Array{Any,1}, ::LinSpace{Float64}) at ./deprecated.jl:644 in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/sparsity_6_l1_recovery/exo4.jl:22 [inlined] in anonymous at ./<missing>:? in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:? in include_string(::String, ::String) at ./loading.jl:441 in execute_request(::ZMQ.Socket, ::IJulia.Msg) at /Users/gpeyre/.julia/v0.5/IJulia/src/execute_request.jl:160 in invokelatest(::Function, ::ZMQ.Socket, ::Vararg{Any,N}) at /Users/gpeyre/.julia/v0.5/Compat/src/Compat.jl:1502 in eventloop(::ZMQ.Socket) at /Users/gpeyre/.julia/v0.5/IJulia/src/eventloop.jl:8 in (::IJulia.##15#21)() at ./task.jl:360 while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/sparsity_6_l1_recovery/exo4.jl, in expression starting on line 11 WARNING: hist(...) and hist!(...) are deprecated. Use fit(Histogram,...) in StatsBase.jl instead. in depwarn(::String, ::Symbol) at ./deprecated.jl:64 in #hist!#994(::Bool, ::Function, ::Array{Int64,1}, ::Array{Any,1}, ::LinSpace{Float64}) at ./deprecated.jl:629 in hist(::Array{Any,1}, ::LinSpace{Float64}) at ./deprecated.jl:644 in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/sparsity_6_l1_recovery/exo4.jl:22 [inlined] in anonymous at ./<missing>:? in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:? in include_string(::String, ::String) at ./loading.jl:441 in execute_request(::ZMQ.Socket, ::IJulia.Msg) at /Users/gpeyre/.julia/v0.5/IJulia/src/execute_request.jl:160 in invokelatest(::Function, ::ZMQ.Socket, ::Vararg{Any,N}) at /Users/gpeyre/.julia/v0.5/Compat/src/Compat.jl:1502 in eventloop(::ZMQ.Socket) at /Users/gpeyre/.julia/v0.5/IJulia/src/eventloop.jl:8 in (::IJulia.##15#21)() at ./task.jl:360 while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions/sparsity_6_l1_recovery/exo4.jl, in expression starting on line 11
## Insert your code here.
Exercise 5
Estimate numerically lower bound on $\de_k^1,\de_k^2$ by Monte-Carlo sampling of sub-matrices.
include("NtSolutions/sparsity_6_l1_recovery/exo5.jl")
## Insert your code here.
We now consider a convolution dictionary $ \Phi $. Such a dictionary is used with sparse regulariz
Second derivative of Gaussian kernel $g$ with a given variance $\si^2$.
sigma = 6
g = x -> (1 - x.^2./sigma^2).*exp(-x.^2./(2*sigma^2))
(::#33) (generic function with 1 method)
Create a matrix $\Phi$ so that $\Phi x = g \star x$ with periodic boundary conditions.
include("ndgrid.jl")
P = 1024
(Y, X) = meshgrid(0:P-1, 0:P-1)
Phi = normalize(g((X - Y + P/2.)%P-P/2.))
could not open file /Users/gpeyre/Dropbox/github/numerical-tours/julia/ndgrid.jl in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
To improve the conditionning of the dictionary, we sub-sample its atoms, so that $ P = \eta N > N $.
eta = 2
N = P/eta
Phi = Phi[:, 1:eta:end]
300×300 Array{Float64,2}: -0.053469 0.0881115 0.0868037 … -0.00996989 -0.0585636 0.0164402 0.0796902 0.0188055 0.025722 -0.0623569 -0.00646314 -0.00191052 -0.0352827 -0.0794237 -0.0370294 -0.017269 0.0359245 0.0216832 -0.0292265 -0.0679947 -0.123559 -0.0629753 0.103768 0.120796 0.0438559 0.0111852 0.0228867 -0.0079232 … 0.0242139 0.0488305 -0.0345084 0.00608496 -0.0685556 -0.0298552 -0.0461946 0.00276266 -0.00457426 0.0113843 0.0393407 0.0983607 -0.0249077 0.0357367 -0.135996 -0.0339714 0.0289777 -0.0181583 -0.0337565 -0.0658213 0.0393399 -0.000954676 -0.0437137 -0.0231178 0.00312412 … 0.0214733 -0.00883125 -0.0419627 -0.0250297 0.0321096 -0.071469 -0.0248835 -0.0845728 -0.037455 -0.0904906 0.0099785 -0.0838841 ⋮ ⋱ 0.143837 -0.0761286 0.0498835 -0.04545 0.0629725 0.0923194 -0.0468798 0.0501927 0.0728482 0.0707002 -0.022314 -0.0734589 0.0360122 … -0.13949 -0.00776711 0.0291248 0.00661235 0.0134569 -0.024864 -0.0150424 -0.0189033 0.0970626 0.0418395 -0.0294127 0.116598 -0.00428868 0.118131 -0.0264479 -0.0932343 -0.154828 0.0778051 0.0158298 0.00234603 0.103604 -0.0754831 -0.101757 0.0495662 -0.021641 … 0.00244892 -0.0197873 -0.0217769 -0.0662395 0.0509761 -0.0209067 0.0325256 -0.0235321 -0.0225997 -0.0209772 -0.0167986 0.0485945 0.0234176 -0.0300984 0.0518727 0.0182376 0.00600304 0.0165291 -0.00807928 -0.0190766 -0.106576 -0.00476109
Plot the correlation function associated to the filter. Can you determine the value of the coherence $\mu(\Phi)$?
c = Phi'*Phi
c = abs(c[:, Int(size(c)[2]/2)])
figure(figsize = (8, 5))
plot(c[Base.div(length(c), 2) - 50 : Base.div(length(c), 2) + 50], ".-")
show()
Create a data a sparse $x_0$ with two diracs of opposite signes with spacing $d$.
twosparse = d -> circshift([1; zeros(d, 1); -1; zeros(Int(N) - d - 2, 1)], round(N/2 - d/2))
(::#35) (generic function with 1 method)
Display $x_0$ and $\Phi x_0$.
x0 = twosparse(50)
figure(figsize = (10, 7))
subplot(2, 1, 1)
plot(x0[:], "r", linewidth = 2)
subplot(2, 1, 2)
plot((Phi*x0)[:], "b", linewidth = 2)
show()
DimensionMismatch("A has dimensions (300,300) but B has dimensions (150,1)") in gemm_wrapper!(::Array{Float64,2}, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at ./linalg/matmul.jl:309 in *(::Array{Float64,2}, ::Array{Float64,2}) at ./linalg/matmul.jl:129
Exercise 6
Plot the evolution of the criteria F, ERC and Coh as a function of $d$. Do the same plot for other signs patterns for $x_0$. Do the same plot for a Dirac comb with a varying spacing $d$.
include("NtSolutions/sparsity_6_l1_recovery/exo6.jl")
# g = (C, I) -> sum(C[:, I], 2)
# figure(figsize = (8, 5))
# dlist = Array{Int64,1}(1 : N/20 - 1)
# criter = zeros(length(dlist), 3)
# for i in 1 : length(dlist)
# s = twosparse(dlist[i])
# I = (supp(s))
# criter[i, :] = [F(Phi, s), erc(Phi,I), werc(Phi,I)]
# end
# criter[criter .< 0] = float("inf")
# plot(dlist, criter, linewidth = 2)
# plot(dlist, dlist.*0 + 1, "k--", linewidth = 2)
# xlim(1, maximum(dlist))
# ylim(minimum(criter), maximum(criter))
# xlabel("d")
# legend(["F", "ERC", "w-ERC"])
# show()
## Insert your code here.