题目：使用l1-magic工具箱求解基追踪(BP)和基追踪降噪(BPDN)

基追踪(Basis Pursuit, BP)和基追踪降噪(Basis PursuitDe-Noising, BPDN)都不能称为是具体的算法，实际上分别是求解一个优化问题：

基追踪：

基追踪降噪：

BP和BPDN可以由多种优化求解方法去解。为了求解这两个优化问题，基追踪可以转化为一个线性规划问题（参见《压缩感知重构算法之基追踪(Basis Pursuit, BP)》，http://blog.csdn.net/jbb0523/article/details/51986554），可以使用MATLAB自带的linprog函数求解；而基追踪降噪可以转化为一个二次规划问题（参见《压缩感知重构算法之基追踪降噪(Basis Pursuit De-Noising, BPDN)》，http://blog.csdn.net/jbb0523/article/details/52013669），可以使用MATLAB自带的quadprog函数求解。

当然也可以用其它方法求解基追踪和基追踪降噪这两个问题。实际上同一个算法，不同的实现方式，运行效率也会有差异。

由Emmanuel Candès and JustinRomberg等人开发的l1-magic工具箱（主页链接：http://users.ece.gatech.edu/justin/l1magic/）可以求解七类优化问题，其中(P1)问题即为基追踪，(P2)问题等价于基追踪降噪，详情可参见工具箱的说明文档^[2]。

(P1)问题：

(P2)问题：

1、(P1)问题

求解(P1)问题的函数是l1eq_pd.m，位于工具箱\l1magic\Optimization目录下。该函数的使用例程参见l1eq_example.m，位于工具箱根目录\l1magic下。在例程l1eq_example中，函数l1eq_pd共有两种调用方式，如下图所示：

默认情况是第一种，第二种为large scale(大规模)，需要注意的是，当使用largescale这种调用方式时，需要额外使用工具箱中的cgsolve.m函数，位于工具箱\l1magic\Optimization目录下。若使用第二种方式，只需将第50行到54行解除注释即可（参见附录代码）。

2、(P2)问题

求解(P2)问题的函数是l1qc_logbarrier.m，位于工具箱\l1magic\Optimization目录下。该函数的使用例程参见l1qc_example.m，位于工具箱根目录\l1magic下。需要注意的是，该函数需要调用l1qc_newton.m函数（位于工具箱\l1magic\Optimization目录下），在例程l1qc_example中，函数l1qc_logbarrier.m共有两种调用方式，如下图所示：

默认情况是第一种，第二种为large scale(大规模)，需要注意的是，当使用largescale这种调用方式时，也需要额外使用工具箱中的cgsolve.m函数。若使用第二种方式，只需将第49行到53行解除注释即可（参见附录代码）。

先简单了解这么一点即可，若实际需要，再详细琢磨一下。

3、总结

总结一下，对于(P1)问题，相关文件有三个：l1eq_pd.m（求解函数）、cgsolve.m（大规模方式时l1eq_pd需调用此文件）、l1eq_example.m（使用范例）；对于(P2)问题，相关文件有四个：l1qc_logbarrier.m（求解函数）、l1qc_newton.m（l1qc_logbarrier需要调用此文件）、cgsolve.m（大规模方式时l1qc_logbarrier需调用此文件）、l1qc_example.m（使用范例）。

值得注意的是，函数l1qc_logbarrier求解基追踪降噪问题的效率和重构效果均要优于本人基于MATLAB自带的quadprog函数编写的BPDN_quadprog（参见《压缩感知重构算法之基追踪降噪(BasisPursuit De-Noising, BPDN)》，http://blog.csdn.net/jbb0523/article/details/52013669），具体原因没有深究。

4、参考文献

【1】ChenS S, Donoho D L, Saunders M A.Atomicdecomposition by basis pursuit[J]. SIAM review, 2001, 43(1): 129-159.(Available at:http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf)

【2】Emmanuel Candès and Justin Romberg.l1-magic : Recovery of Sparse Signals via Convex Programming[EB/OL].http://users.ece.gatech.edu/justin/l1magic/downloads/l1magic.pdf

5、附录

作为一种备份，也作为一篇文档的完整性，这里将(P1)(P2)所涉及的MATLAB文件附在此，如果不愿意自己去下载，可以将这些代码保存为MATLAB文件即可。

(1)第一个是两个问题在large scale方式下都要调用的cgsolve.m

% cgsolve.m
%
% Solve a symmetric positive definite system Ax = b via conjugate gradients.
%
% Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
%
% A - Either an NxN matrix, or a function handle.
%
% b - N vector
%
% tol - Desired precision.  Algorithm terminates when
%    norm(Ax-b)/norm(b) < tol .
%
% maxiter - Maximum number of iterations.
%
% verbose - If 0, do not print out progress messages.
%    If and integer greater than 0, print out progress every 'verbose' iters.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)if (nargin < 5), verbose = 1; endimplicit = isa(A,'function_handle');x = zeros(length(b),1);
r = b;
d = r;
delta = r'*r;
delta0 = b'*b;
numiter = 0;
bestx = x;
bestres = sqrt(delta/delta0);
while ((numiter < maxiter) && (delta > tol^2*delta0))% q = A*dif (implicit), q = A(d);  else  q = A*d;  endalpha = delta/(d'*q);x = x + alpha*d;if (mod(numiter+1,50) == 0)% r = b - Aux*xif (implicit), r = b - A(x);  else  r = b - A*x;  endelser = r - alpha*q;enddeltaold = delta;delta = r'*r;beta = delta/deltaold;d = r + beta*d;numiter = numiter + 1;if (sqrt(delta/delta0) < bestres)bestx = x;bestres = sqrt(delta/delta0);end    if ((verbose) && (mod(numiter,verbose)==0))disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ...numiter, bestres, sqrt(delta/delta0)));endendif (verbose)disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres));
end
x = bestx;
res = bestres;
iter = numiter;

(2)第二个是(P1)问题求解函数l1eq_pd.m

% l1eq_pd.m
%
% Solve
% min_x ||x||_1  s.t.  Ax = b
%
% Recast as linear program
% min_{x,u} sum(u)  s.t.  -u <= x <= u,  Ax=b
% and use primal-dual interior point method
%
% Usage: xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)
%
% x0 - Nx1 vector, initial point.
%
% A - Either a handle to a function that takes a N vector and returns a K
%     vector , or a KxN matrix.  If A is a function handle, the algorithm
%     operates in "largescale" mode, solving the Newton systems via the
%     Conjugate Gradients algorithm.
%
% At - Handle to a function that takes a K vector and returns an N vector.
%      If A is a KxN matrix, At is ignored.
%
% b - Kx1 vector of observations.
%
% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if
%     the duality gap is less than pdtol).
%     Default = 1e-3.
%
% pdmaxiter - Maximum number of primal-dual iterations.
%     Default = 50.
%
% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
%     Default = 1e-8.
%
% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
%     if A is a matrix.
%     Default = 200.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%function xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)largescale = isa(A,'function_handle');if (nargin < 5), pdtol = 1e-3;  end
if (nargin < 6), pdmaxiter = 50;  end
if (nargin < 7), cgtol = 1e-8;  end
if (nargin < 8), cgmaxiter = 200;  endN = length(x0);alpha = 0.01;
beta = 0.5;
mu = 10;gradf0 = [zeros(N,1); ones(N,1)];% starting point --- make sure that it is feasible
if (largescale)if (norm(A(x0)-b)/norm(b) > cgtol)disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');AAt = @(z) A(At(z));[w, cgres, cgiter] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);if (cgres > 1/2)disp('A*At is ill-conditioned: cannot find starting point');xp = x0;return;endx0 = At(w);end
elseif (norm(A*x0-b)/norm(b) > cgtol)disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');opts.POSDEF = true; opts.SYM = true;[w, hcond] = linsolve(A*A', b, opts);if (hcond < 1e-14)disp('A*At is ill-conditioned: cannot find starting point');xp = x0;return;endx0 = A'*w;end
end
x = x0;
u = (0.95)*abs(x0) + (0.10)*max(abs(x0));% set up for the first iteration
fu1 = x - u;
fu2 = -x - u;
lamu1 = -1./fu1;
lamu2 = -1./fu2;
if (largescale)v = -A(lamu1-lamu2);Atv = At(v);rpri = A(x) - b;
elsev = -A*(lamu1-lamu2);Atv = A'*v;rpri = A*x - b;
endsdg = -(fu1'*lamu1 + fu2'*lamu2);
tau = mu*2*N/sdg;rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
resnorm = norm([rdual; rcent; rpri]);pditer = 0;
done = (sdg < pdtol) | (pditer >= pdmaxiter);
while (~done)pditer = pditer + 1;w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv;w2 = -1 - 1/tau*(1./fu1 + 1./fu2);w3 = -rpri;sig1 = -lamu1./fu1 - lamu2./fu2;sig2 = lamu1./fu1 - lamu2./fu2;sigx = sig1 - sig2.^2./sig1;if (largescale)w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1));h11pfun = @(z) -A(1./sigx.*At(z));[dv, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);if (cgres > 1/2)disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');xp = x;returnenddx = (w1 - w2.*sig2./sig1 - At(dv))./sigx;Adx = A(dx);Atdv = At(dv);elsew1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1)));H11p = A*(sparse(diag(1./sigx))*A');opts.POSDEF = true; opts.SYM = true;[dv,hcond] = linsolve(H11p, w1p, opts);if (hcond < 1e-14)disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');xp = x;returnenddx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx;Adx = A*dx;Atdv = A'*dv;enddu = (w2 - sig2.*dx)./sig1;dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1;dlamu2 = (lamu2./fu2).*(dx+du) - lamu2 - 1/tau*1./fu2;% make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0indp = find(dlamu1 < 0);  indn = find(dlamu2 < 0);s = min([1; -lamu1(indp)./dlamu1(indp); -lamu2(indn)./dlamu2(indn)]);indp = find((dx-du) > 0);  indn = find((-dx-du) > 0);s = (0.99)*min([s; -fu1(indp)./(dx(indp)-du(indp)); -fu2(indn)./(-dx(indn)-du(indn))]);% backtracking line searchsuffdec = 0;backiter = 0;while (~suffdec)xp = x + s*dx;  up = u + s*du; vp = v + s*dv;  Atvp = Atv + s*Atdv; lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;fu1p = xp - up;  fu2p = -xp - up;  rdp = gradf0 + [lamu1p-lamu2p; -lamu1p-lamu2p] + [Atvp; zeros(N,1)];rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau);rpp = rpri + s*Adx;suffdec = (norm([rdp; rcp; rpp]) <= (1-alpha*s)*resnorm);s = beta*s;backiter = backiter + 1;if (backiter > 32)disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')xp = x;returnendend% next iterationx = xp;  u = up;v = vp;  Atv = Atvp; lamu1 = lamu1p;  lamu2 = lamu2p;fu1 = fu1p;  fu2 = fu2p;% surrogate duality gapsdg = -(fu1'*lamu1 + fu2'*lamu2);tau = mu*2*N/sdg;rpri = rpp;rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];resnorm = norm([rdual; rcent; rpri]);done = (sdg < pdtol) | (pditer >= pdmaxiter);disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e, Primal res = %8.3e',...pditer, tau, sum(u), sdg, norm(rdual), norm(rpri)));if (largescale)disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));elsedisp(sprintf('                  H11p condition number = %8.3e', hcond));endend

(3)第三个是 (P1)问题使用范例l1eq_example.m（处于large scale方式下）

% l1eq_example.m
%
% Test out l1eq code (l1 minimization with equality constraints).
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%% put key subdirectories in path if not already there
% path(path, './Optimization');
% path(path, './Data');% To reproduce the example in the documentation, uncomment the
% two lines below
%load RandomStates
%rand('state', rand_state);
%randn('state', randn_state);% signal length
N = 512;
% number of spikes in the signal
T = 20;
% number of observations to make
K = 120;% random +/- 1 signal
x = zeros(N,1);
q = randperm(N);
x(q(1:T)) = sign(randn(T,1));% measurement matrix
disp('Creating measurment matrix...');
A = randn(K,N);
A = orth(A')';
disp('Done.');% observations
y = A*x;% initial guess = min energy
x0 = A'*y;% solve the LP
% tic
% xp = l1eq_pd(x0, A, [], y, 1e-3);
% toc% large scale
Afun = @(z) A*z;
Atfun = @(z) A'*z;
tic
xp = l1eq_pd(x0, Afun, Atfun, y, 1e-3, 30, 1e-8, 200);
toc

(4)第四个是 (P2)问题求解函数l1qc_logbarrier.m

% l1qc_logbarrier.m
%
% Solve quadratically constrained l1 minimization:
% min ||x||_1   s.t.  ||Ax - b||_2 <= \epsilon
%
% Reformulate as the second-order cone program
% min_{x,u}  sum(u)   s.t.    x - u <= 0,
%                            -x - u <= 0,
%      1/2(||Ax-b||^2 - \epsilon^2) <= 0
% and use a log barrier algorithm.
%
% Usage:  xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
%
% x0 - Nx1 vector, initial point.
%
% A - Either a handle to a function that takes a N vector and returns a K
%     vector , or a KxN matrix.  If A is a function handle, the algorithm
%     operates in "largescale" mode, solving the Newton systems via the
%     Conjugate Gradients algorithm.
%
% At - Handle to a function that takes a K vector and returns an N vector.
%      If A is a KxN matrix, At is ignored.
%
% b - Kx1 vector of observations.
%
% epsilon - scalar, constraint relaxation parameter
%
% lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
%         Also, the number of log barrier iterations is completely
%         determined by lbtol.
%         Default = 1e-3.
%
% mu - Factor by which to increase the barrier constant at each iteration.
%      Default = 10.
%
% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
%     Default = 1e-8.
%
% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
%     if A is a matrix.
%     Default = 200.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)  largescale = isa(A,'function_handle');if (nargin < 6), lbtol = 1e-3; end
if (nargin < 7), mu = 10; end
if (nargin < 8), cgtol = 1e-8; end
if (nargin < 9), cgmaxiter = 200; endnewtontol = lbtol;
newtonmaxiter = 50;N = length(x0);% starting point --- make sure that it is feasible
if (largescale)if (norm(A(x0)-b) > epsilon)disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');AAt = @(z) A(At(z));[w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);if (cgres > 1/2)disp('A*At is ill-conditioned: cannot find starting point');xp = x0;return;endx0 = At(w);end
elseif (norm(A*x0-b) > epsilon)disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');opts.POSDEF = true; opts.SYM = true;[w, hcond] = linsolve(A*A', b, opts);if (hcond < 1e-14)disp('A*At is ill-conditioned: cannot find starting point');xp = x0;return;endx0 = A'*w;end
end
x = x0;
u = (0.95)*abs(x0) + (0.10)*max(abs(x0));disp(sprintf('Original l1 norm = %.3f, original functional = %.3f', sum(abs(x0)), sum(u)));% choose initial value of tau so that the duality gap after the first
% step will be about the origial norm
tau = max((2*N+1)/sum(abs(x0)), 1);lbiter = ceil((log(2*N+1)-log(lbtol)-log(tau))/log(mu));
disp(sprintf('Number of log barrier iterations = %d\n', lbiter));totaliter = 0;for ii = 1:lbiter[xp, up, ntiter] = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter);totaliter = totaliter + ntiter;disp(sprintf('\nLog barrier iter = %d, l1 = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...ii, sum(abs(xp)), sum(up), tau, totaliter));x = xp;u = up;tau = mu*tau;end

(5)第五个是 (P2)问题求解函数须调用的l1qc_newton.m

% l1qc_newton.m
%
% Newton algorithm for log-barrier subproblems for l1 minimization
% with quadratic constraints.
%
% Usage:
% [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau,
%                             newtontol, newtonmaxiter, cgtol, cgmaxiter)
%
% x0,u0 - starting points
%
% A - Either a handle to a function that takes a N vector and returns a K
%     vector , or a KxN matrix.  If A is a function handle, the algorithm
%     operates in "largescale" mode, solving the Newton systems via the
%     Conjugate Gradients algorithm.
%
% At - Handle to a function that takes a K vector and returns an N vector.
%      If A is a KxN matrix, At is ignored.
%
% b - Kx1 vector of observations.
%
% epsilon - scalar, constraint relaxation parameter
%
% tau - Log barrier parameter.
%
% newtontol - Terminate when the Newton decrement is <= newtontol.
%         Default = 1e-3.
%
% newtonmaxiter - Maximum number of iterations.
%         Default = 50.
%
% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
%     Default = 1e-8.
%
% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
%     if A is a matrix.
%     Default = 200.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%function [xp, up, niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) % check if the matrix A is implicit or explicit
largescale = isa(A,'function_handle');% line search parameters
alpha = 0.01;
beta = 0.5;  if (~largescale), AtA = A'*A; end% initial point
x = x0;
u = u0;
if (largescale), r = A(x) - b; else  r = A*x - b; end
fu1 = x - u;
fu2 = -x - u;
fe = 1/2*(r'*r - epsilon^2);
f = sum(u) - (1/tau)*(sum(log(-fu1)) + sum(log(-fu2)) + log(-fe));niter = 0;
done = 0;
while (~done)if (largescale), atr = At(r); else  atr = A'*r; endntgz = 1./fu1 - 1./fu2 + 1/fe*atr;ntgu = -tau - 1./fu1 - 1./fu2;gradf = -(1/tau)*[ntgz; ntgu];sig11 = 1./fu1.^2 + 1./fu2.^2;sig12 = -1./fu1.^2 + 1./fu2.^2;sigx = sig11 - sig12.^2./sig11;w1p = ntgz - sig12./sig11.*ntgu;if (largescale)h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr;[dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);if (cgres > 1/2)disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');xp = x;  up = u;returnendAdx = A(dx);elseH11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr';opts.POSDEF = true; opts.SYM = true;[dx,hcond] = linsolve(H11p, w1p, opts);if (hcond < 1e-14)disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');xp = x;  up = u;returnendAdx = A*dx;enddu = (1./sig11).*ntgu - (sig12./sig11).*dx;  % minimum step size that stays in the interiorifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0);aqe = Adx'*Adx;   bqe = 2*r'*Adx;   cqe = r'*r - epsilon^2;smax = min(1,min([...-fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ...(-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe)]));s = (0.99)*smax;% backtracking line searchsuffdec = 0;backiter = 0;while (~suffdec)xp = x + s*dx;  up = u + s*du;  rp = r + s*Adx;fu1p = xp - up;  fu2p = -xp - up;  fep = 1/2*(rp'*rp - epsilon^2);fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep));flin = f + alpha*s*(gradf'*[dx; du]);suffdec = (fp <= flin);s = beta*s;backiter = backiter + 1;if (backiter > 32)disp('Stuck on backtracking line search, returning previous iterate.  (See Section 4 of notes for more information.)');xp = x;  up = u;returnendend% set up for next iterationx = xp; u = up;  r = rp;fu1 = fu1p;  fu2 = fu2p;  fe = fep;  f = fp;lambda2 = -(gradf'*[dx; du]);stepsize = s*norm([dx; du]);niter = niter + 1;done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ...niter, f, lambda2/2, stepsize));if (largescale)disp(sprintf('                CG Res = %8.3e, CG Iter = %d', cgres, cgiter));elsedisp(sprintf('                  H11p condition number = %8.3e', hcond));endend

(6)第六个是(P2)问题使用范例l1qc_example.m（处于large scale方式下）

% l1qc_example.m
%
% Test out l1qc code (l1 minimization with quadratic constraint).
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%% put optimization code in path if not already there
path(path, './Optimization');% signal length
N = 512;% number of spikes to put down
T = 20;% number of observations to make
K = 120;% random +/- 1 signal
x = zeros(N,1);
q = randperm(N);
x(q(1:T)) = sign(randn(T,1));% measurement matrix
disp('Creating measurment matrix...');
A = randn(K,N);
A = orth(A')';
disp('Done.');% noisy observations
sigma = 0.005;
e = sigma*randn(K,1);
y = A*x + e;% initial guess = min energy
x0 = A'*y;% take epsilon a little bigger than sigma*sqrt(K)
epsilon =  sigma*sqrt(K)*sqrt(1 + 2*sqrt(2)/sqrt(K));% tic
% xp = l1qc_logbarrier(x0, A, [], y, epsilon, 1e-3);
% toc% large scale
Afun = @(z) A*z;
Atfun = @(z) A'*z;
tic
xp = l1qc_logbarrier(x0, Afun, Atfun, y, epsilon, 1e-3, 50, 1e-8, 500);
toc

使用l1-magic工具箱求解基追踪(BP)和基追踪降噪(BPDN)相关推荐

matlab中PDE工具箱如何使用,使用PDE工具箱求解偏微分方程
在科学技术各领域中,有很多问题都可以归结为偏微分方程问题.在物理专业的力学.热学.电学.光学.近代物理课程中都可遇见偏微分方程. 偏微分方程,再加上边界条件.初始条件构成的数学模型,只有在很特殊情况下 ...
matlab偏微分方程工具箱求解
Matlab的偏微分方程工具箱求解方法这一节我们主要用matlab自带的偏微分方程的工具箱函数求解一.偏微分方程组的matlab求解语句该命令用以求解以下的PDEPDEPDE方程式: c(x ...
Magic Leap开发指南（8）-- 眼球追踪(Lumin Runtime)
上一篇教程(Magic Leap开发指南(7)-- 眼球追踪(Unity))我们了解了如何在Unity中使用Eye Tracking来完成一些小项目,这篇我们继续通过Lumin Runtime来运用眼 ...
【运筹学】人工变量法总结 ( 人工变量法解的分析 | 标准型变换 | 构造单位阵 | 目标函数引入 M | 计算检验数 | 选择入基变量 | 选择出基变量 | 中心元变换 | ) ★★
文章目录一.人工变量法及解的分析二.案例三.线性规划标准型变换四.人工变量法构造单位阵五.初始单纯形表六.初始单纯形表 : 计算非基变量检验数七.初始单纯形表 : 最优解判定八.初始单 ...
【运筹学】线性规划单纯形法案例二 ( 案例解析 | 标准形转化 | 查找初始基可行解 | 最优解判定 | 查找入基变量与出基变量 | 第一次迭代 )
文章目录一.线性规划示例二.转化成标准形式三.初始基可行解四.列出单纯形表五.计算检验数六.选择入基变量与出基变量七.第一次迭代 : 列出单纯形表一.线性规划示例线性规划示例 : 使 ...
【运筹学】线性规划数学模型 ( 单纯形法 | 第一次迭代 | 方程组同解变换 | 计算新单纯形表 | 计算检验数 | 入基变量选择 | 出基变量选择 )
文章目录一.初始基可行解后第一次迭代二.迭代后新的单纯形表三.方程组同解变换四.生成新的单纯形表五.解出基可行解六.计算检验数 σj\sigma_jσj 并选择入基变量七.计算 θ\t ...
【运筹学】线性规划人工变量法 ( 人工变量法案例 | 第二次迭代 | 中心元变换 | 检验数计算 | 最优解判定 | 选择入基变量 | 选择出基变量 )
文章目录一.第二次迭代 : 中心元变换二.第二次迭代 : 单纯形表三.第二次迭代 : 计算检验数四.第二次迭代 : 最优解判定五.第二次迭代 : 选择入基变量六.第二次迭代 : 选择出基变 ...
【运筹学】线性规划人工变量法 ( 人工变量法案例 | 第一次迭代 | 中心元变换 | 检验数计算 | 选择入基变量 | 选择出基变量 )
文章目录一.第一次迭代 : 中心元变换二.第一次迭代 : 单纯形表三.第一次迭代 : 计算检验数四.第一次迭代 : 最优解判定五.第一次迭代 : 选择入基变量六.第一次迭代 : 选择出基变 ...
unity怎么实现人脸追踪_Unity 2019.2 beta为AR增加面部追踪、2D图像追踪、3D对象追踪等功能...
Unity今天正式放出了Unity 2019.2 beta.对于这个版本,Unity集成了热门的Polybrush工具,添加了Unity Distribution Platform,同时扩展了用于XR ...

使用l1-magic工具箱求解基追踪(BP)和基追踪降噪(BPDN)