📄 index.html
字号:
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"><html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Sparse Spikes Deconvolution with L1 Pursuit</title> <meta name="generator" content="MATLAB 7.4"> <meta name="date" content="2008-10-15"> <meta name="m-file" content="index"> <LINK REL="stylesheet" HREF="style.css" TYPE="text/css"> </head> <body> <div class="content"> <h1>Sparse Spikes Deconvolution with L1 Pursuit</h1> <introduction> <p>This numerical tour explores the use of L1 regularization (basis pursuit denoising) to solve seismic sparse spikes deconvolution.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#3">Installing toolboxes and setting up the path.</a></li> <li><a href="#10">Seismic Wavelets and Seismic Imaging</a></li> <li><a href="#16">Homotopy Method for L1 Regularization</a></li> <li><a href="#36">Iterative Soft Thresholding for L1 Minimization</a></li> </ul> </div> <p>Sparse spikes deconvolution is one of the oldest inverse problems, that is a stilized version of recovery in seismic imaging. The ground is modeled as a 1D profile <tt>x</tt>, mostly zeros with a few spikes accounting for interfaces between layers in the ground. The observation is the convolution <tt>D*x</tt> of <tt>x</tt> against a wavelet filter, where <tt>D(:,1)</tt> is the basis wavelet function (<tt>D</tt> is the convolution operator. </p> <p>The goal of sparse spike deconvolution is to recover an approximation of <tt>x</tt> given noisy measurement <tt>y=D*x+w</tt>. Since the convolution destroys many low and high frequencies, this requires some prior information to regularize the inverse problem. Since <tt>x</tt> is composed of a few spikes, the sparsity is a wonderful prior, that is exploited either by L1 minimization (basis pursuit) or by non-linear greedy processes (matching pursuits). </p> <h2>Installing toolboxes and setting up the path.<a name="3"></a></h2> <p>You need to download the <a href="../toolbox_general.zip">general purpose toolbox</a> and the <a href="../toolbox_signal.zip">signal toolbox</a>. </p> <p>You need to unzip these toolboxes in your working directory, so that you have <tt>toolbox_general/</tt> and <tt>toolbox_signal/</tt> in your directory. </p> <p><b>For Scilab user:</b> you must replace the Matlab comment '%' by its Scilab counterpart '//'. </p> <p><b>Recommandation:</b> You should create a text file named for instance <tt>numericaltour.sce</tt> (in Scilabe) or <tt>numericaltour.m</tt> to write all the Scilab/Matlab command you want to execute. Then, simply run <tt>exec('numericaltour.sce');</tt> (in Scilab) or <tt>numericaltour;</tt> (in Matlab) to run the commands. </p> <p>Execute this line only if you are using Matlab.</p><pre class="codeinput">getd = @(p)path(path,p); <span class="comment">% scilab users must *not* execute this</span></pre><p>Then you can add these toolboxes to the path.</p><pre class="codeinput"><span class="comment">% Add some directories to the path</span>getd(<span class="string">'toolbox_signal/'</span>);getd(<span class="string">'toolbox_general/'</span>);</pre><h2>Seismic Wavelets and Seismic Imaging<a name="10"></a></h2> <p>We first build the dictionary <tt>D</tt> used for seismic imaging. </p><pre class="codeinput">set_rand_seeds(123456,21);</pre><p>First we load a seismic filter, which is a second derivative of a Gaussian.</p><pre class="codeinput"><span class="comment">% dimension of the signal</span>n = 1024;<span class="comment">% width of the filter</span>s = 13;<span class="comment">% second derivative of Gaussian</span>t = -n/2:n/2-1;h = (1-t.^2/s^2).*exp( -(t.^2)/(2*s^2) );h = h-mean(h);<span class="comment">% normalize it</span>h = h/norm(h);<span class="comment">% recenter the filter for periodic boundary conditions</span>h1 = fftshift(h);</pre><p>We compute the filtering matrix. To stabilize the recovery, we sub-sample by a factor of 2 the filtering.</p><pre class="codeinput"><span class="comment">% sub-sampling (distance between wavelets)</span>sub = 2;<span class="comment">% number of atoms in the dictionary</span>p = n/sub;<span class="comment">% the dictionary, with periodic boundary conditions</span>[Y,X] = meshgrid(1:sub:n,1:n);D = reshape( h1(mod(X-Y,n)+1), [n p]);</pre><p>We generate a sparse signal to recover.</p><pre class="codeinput"><span class="comment">% spacing min and max between the spikes.</span>m = 5; M = 40;k = floor( (p+M)*2/(M+m) )-1;spc = linspace(M,m,k)';<span class="comment">% location of the spikes</span>sel = round( cumsum(spc) );sel(sel>p) = [];<span class="comment">% randomization of the signs and values</span>x = zeros(p,1);si = (-1).^(1:length(sel))'; si = si(randperm(length(si)));<span class="comment">% creating of the sparse spikes signal.</span>x(sel) = si;x = x .* (1-rand(p,1)*.5);<span class="comment">% sparsity of the solution</span>M = sum(x~=0);</pre><p>Now we perform the measurements.</p><pre class="codeinput"><span class="comment">% noise level</span>sigma = .06*max(h);<span class="comment">% noise</span>w = randn(n,1)*sigma;<span class="comment">% measures</span>y = D*x + w;</pre><p>Display.</p><pre class="codeinput">clf;subplot(2,1,1);plot_sparse_diracs(x);title(<span class="string">'Signal x'</span>);subplot(2,1,2);plot(y);set_graphic_sizes([], 20);axis(<span class="string">'tight'</span>);title(<span class="string">'Noisy measurements y=D*x+w'</span>);</pre><img vspace="5" hspace="5" src="index_01.png"> <h2>Homotopy Method for L1 Regularization<a name="16"></a></h2> <p>Instead of using a greedy matching pursuit to recover an approximation of <tt>x</tt>, we now turn to more global minimization method. L1 minimization is a convex regularization that approximates the L0 pseudo-norm regularization (that would be ideal to find highly sparse Diracs train). The minimization reads </p> <p><tt>x(lambda) = argmin_x 1/2*norm(y-D*x)^2 + lambda*norm(x,1)</tt></p> <p>The Homotopy method starts with a large <tt>lambda</tt> so that <tt>x(lambda)=0</tt> and progressively decays the value of <tt>lambda</tt> so that the sparsity of <tt>x(lambda)</tt> only increases/decreases by 1 at each step. </p> <p>The intial solution is 0.</p><pre class="codeinput">xbp = zeros(p,1);</pre><p>The solution <tt>xbp</tt> is update <tt>xbp = xbp + gamma*d</tt> where the direction of update is supported on the set of vectors that maximally correlate with the residual <tt>y-D*xbp</tt></p><pre class="codeinput"><span class="comment">% compute the correlation with the residual</span>C = D'*(y-D*xbp);lambda = max(abs(C));<span class="comment">% find the locations that maximally correlates</span>S = find( abs(abs(C/lambda)-1)<1e-9);<span class="comment">% compute the complementary set</span>I = ones(p,1); I(S)=0;Sc = find(I);</pre><p>The direction of descent <tt>d(S)</tt> on <tt>S</tt> is computed so that its image by <tt>D(:,S)*d(S)</tt> correlates as +1 or -1 with the atoms in <tt>S</tt> (same speed on all the coefficients). It is zero outside <tt>S</tt>. </p><pre class="codeinput">d = zeros(p,1);d(S) = (D(:,S)'*D(:,S)) \ sign( C(S) );</pre><p>Now we compute the value <tt>gamma</tt> so that either </p> <div> <ul> <li>1) <tt>xbp + gamma*d</tt> correlates as much as <tt>lambda</tt> with one atom outside <tt>S</tt>. </li> <li>2) <tt>xbp + gamma*d</tt> correlates as much as <tt>-lambda</tt> with one atom outside <tt>S</tt>. </li> <li>3) one of the coordinates <tt>xbp + gamma*d</tt> in <tt>S</tt> becomes 0. </li> </ul> </div> <p>In situations 1) and 2), the number of non-zero coordinates of <tt>xbp</tt> (its sparsity) increases by 1. In situation 3), its sparsity stays constant because one coefficients appears and another deasapears. </p> <p>Compute minimum <tt>gamma</tt> so that situation 1) is in force. </p><pre class="codeinput">v = D(:,S)*d(S);w = ( lambda-C(Sc) ) ./ ( 1 - D(:,Sc)'*v );gamma1 = min(w(w>0));</pre><p>Compute minimum <tt>gamma</tt> so that situation 2) is in force. </p><pre class="codeinput">w = ( lambda+C(Sc) ) ./ ( 1 + D(:,Sc)'*v );gamma2 = min(w(w>0));</pre><p>Compute minimum <tt>gamma</tt> so that situation 3) is in force. </p><pre class="codeinput">w = -xbp(S)./d(S);gamma3 = min(w(w>0));</pre><p>Compute minimum <tt>gamma</tt> so that 1), 2) or 3) is in force, and update the solution. </p><pre class="codeinput">gamma = min([gamma1 gamma2 gamma3]);</pre><p>We can display the residual correlation before/after the update. You see that a new maximally correlating atoms has appeared, and that the overall correlation has decayed a little. </p><pre class="codeinput">clf;subplot(2,1,1);plot( abs(D'*(y-D*xbp)), <span class="string">'.-'</span> ); axis <span class="string">tight</span>;title(<span class="string">'Correlation before update'</span>);subplot(2,1,2);plot( abs(D'*(y-D*(xbp+gamma*d))), <span class="string">'.-'</span> ); axis <span class="string">tight</span>;title(<span class="string">'Correlation after update'</span>);</pre><img vspace="5" hspace="5" src="index_02.png"> <p>Update the solution.</p><pre class="codeinput">xbp = xbp + gamma*d;</pre><p><i>Exercice 1:</i> (the solution is <a href="../private/sparsity_seismic_bp/exo1.m">exo1.m</a>) Implements the homotopy algorithm by applying several times (let's say up to <tt>1.5*M</tt> times) the previous update rule. Keep track of the evolution of <tt>lambda</tt> and of the evolution of the sparsity. </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_03.png"> <p>Display the solution computed by L1.</p><pre class="codeinput">clf;subplot(2,1,1);plot_sparse_diracs(x);title(<span class="string">'Signal x'</span>);subplot(2,1,2);plot_sparse_diracs(xbp);title(<span class="string">'Recovered by L1 minimization'</span>);</pre><img vspace="5" hspace="5" src="index_04.png"> <p>The solution computed by L1 as a tendency to under-estimate the true values of the coefficients. This is bias inherent to L1 minimization. To remove this bias, one can perform a back projection that performs an L2 best fit on the support selected by L1. </p><pre class="codeinput"><span class="comment">% find the support</span>sel = find(xbp~=0);<span class="comment">% perform the fit</span>xproj = zeros(p,1);xproj(sel) = D(:,sel) \ y;<span class="comment">% display</span>clf;subplot(2,1,1);plot_sparse_diracs(xbp);title(<span class="string">'Recovered by L1 minimization'</span>);subplot(2,1,2);plot_sparse_diracs(xproj);title(<span class="string">'Back-projection'</span>);</pre><img vspace="5" hspace="5" src="index_05.png"> <p><i>Exercice 2:</i> (the solution is <a href="../private/sparsity_seismic_bp/exo2.m">exo2.m</a>) In the same code, add a tracking of the error <tt>norm(x-xproj)</tt>, and return the solution that mimizes this error. Save the value of the optimal <tt>lambda</tt> in <tt>lambda_opt</tt> and the value <tt>xbp</tt> of the solution for this <tt>lambda</tt> for a future use. </p><pre class="codeinput">exo2;</pre><img vspace="5" hspace="5" src="index_06.png"> <p>Display the best solution</p><pre class="codeinput">err_bp = norm(x-xproj)/norm(x);clf;subplot(2,1,1);plot_sparse_diracs(x);title(<span class="string">'Signal x'</span>);subplot(2,1,2);plot_sparse_diracs(xproj);title( strcat([<span class="string">'L1 recovery, error = '</span> num2str(err_bp, 3)]));</pre><img vspace="5" hspace="5" src="index_07.png"> <h2>Iterative Soft Thresholding for L1 Minimization<a name="36"></a></h2> <p>If the value of <tt>lambda</tt> is known, one can use an iterative algorithm to find the solution <tt>xbp</tt> of the L1 minimization. </p> <p>We use the optimal value of <tt>lambda</tt> already computed. We also save the true solution computed by homotopy. </p><pre class="codeinput">lambda = lambda_opt;xbp_opt = xbp;</pre><p>The gradient descent step size depends on the conditionning of the matrix.</p><pre class="codeinput">tau = 2/norm(D).^2;</pre><p>The iterative algorithm starts with the zero vector.</p><pre class="codeinput">xbp = zeros(p,1);</pre><p>The gradient step updates the value of the solution by decaying the value of <tt>norm(y-D*xbp)^2</tt>. </p><pre class="codeinput">xbp = xbp + tau*D'*( y-D*xbp );</pre><p>The thresholding step improves the sparsity of the solution.</p><pre class="codeinput">xbp_prev = xbp;xbp = perform_thresholding( xbp, tau*lambda, <span class="string">'soft'</span> );<span class="comment">% display</span>clf;subplot(2,1,1);plot(xbp_prev); axis(<span class="string">'tight'</span>);set_graphic_sizes([], 20);title(<span class="string">'Estimate before soft thresholding'</span>);subplot(2,1,2);plot(xbp); axis(<span class="string">'tight'</span>);set_graphic_sizes([], 20);title(<span class="string">'Estimate after soft thresholding'</span>);</pre><img vspace="5" hspace="5" src="index_08.png"> <p><i>Exercice 3:</i> (the solution is <a href="../private/sparsity_seismic_bp/exo3.m">exo3.m</a>) Implement the iterative soft thresholding algorithm by applying many times these two steps. Keep track of the variational energy minimized by this algorithm. Keep also track of the distance to the solution. </p><pre class="codeinput">exo3;</pre><img vspace="5" hspace="5" src="index_09.png"> <p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Sparse Spikes Deconvolution with L1 Pursuit% This numerical tour explores the use of L1 regularization (basis pursuit denoising) to solve% seismic sparse spikes deconvolution.%%% Sparse spikes deconvolution is one of the oldest inverse problems, that% is a stilized version of recovery in seismic imaging.% The ground is modeled as a 1D profile |x|, mostly zeros with a few spikes% accounting for interfaces between layers in the ground. % The observation is the convolution |D*x| of |x| against a wavelet filter,% where |D(:,1)| is the basis wavelet function (|D| is the convolution% operator. %%% The goal of sparse spike deconvolution is to recover an approximation% of |x| given noisy measurement |y=D*x+w|. Since the convolution destroys% many low and high frequencies, this requires some prior information to% regularize the inverse problem. Since |x| is composed of a few spikes,% the sparsity is a wonderful prior, that is exploited either by L1% minimization (basis pursuit) or by non-linear greedy processes (matching% pursuits).%% Installing toolboxes and setting up the path.%%% You need to download the % <../toolbox_general.zip general purpose toolbox>% and the <../toolbox_signal.zip signal toolbox>.%%% You need to unzip these toolboxes in your working directory, so% that you have |toolbox_general/| and |toolbox_signal/| in your directory.%%% *For Scilab user:* you must replace the Matlab comment '%' by its Scilab% counterpart '//'.%%% *Recommandation:* You should create a text file named for instance% |numericaltour.sce| (in Scilabe) or |numericaltour.m| to write all the% Scilab/Matlab command you want to execute. Then, simply run% |exec('numericaltour.sce');| (in Scilab) or |numericaltour;| (in Matlab)% to run the commands. %%% Execute this line only if you are using Matlab.getd = @(p)path(path,p); % scilab users must *not* execute this%%% Then you can add these toolboxes to the path.% Add some directories to the pathgetd('toolbox_signal/');getd('toolbox_general/');%% Seismic Wavelets and Seismic Imaging% We first build the dictionary |D| used for seismic imaging.⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -