bpdq_proj_lpball : Projection onto Lp ball via Newton's method function [x_out,n] = bpdq_proj_lpball(y_in,c,r,p) This function computes the projection of y_in onto the lp ball, ie solves x_out = argmin ||x-y_in|| s.t. ||x-c||_p <= r Algorithm proceeds by multivariate Newton's method applied to find the zeros of the N+1 dimensional nonlinear equation F(z)=0 arising from lagrange multiplier condition for the constrained optimization problem. Algorithm should not be used for p<2 as Jacobian entries contain z_k^(p-2) ... Inputs : y_in - input vector, to be projected c - center of Lp ball r - radius of Lp ball p - exponent for Lp ball Outputs: x_out - computed projection n - number of iterations This is a matlab implementation, if possible the compiled c++ implementation bpdq_proj_lpball_mex should be used instead. This matlab file may be considered a more easily readable reference implementation of the algorithm, and is included largely for pedagogic purposes This file is part of BPDQ Toolbox (Basis Pursuit DeQuantizer) Copyright (C) 2009, the BPDQ Team (see the file AUTHORS distributed with this library) (See the notice at the end of the file.)
0001 % bpdq_proj_lpball : Projection onto Lp ball via Newton's method 0002 % 0003 % function [x_out,n] = bpdq_proj_lpball(y_in,c,r,p) 0004 % 0005 % This function computes the projection of y_in onto the lp ball, ie 0006 % solves 0007 % 0008 % x_out = argmin ||x-y_in|| s.t. ||x-c||_p <= r 0009 % 0010 % Algorithm proceeds by multivariate Newton's method applied to 0011 % find the zeros of the N+1 dimensional nonlinear equation F(z)=0 0012 % arising from lagrange multiplier condition for the constrained 0013 % optimization problem. 0014 % 0015 % Algorithm should not be used for p<2 as Jacobian entries contain 0016 % z_k^(p-2) ... 0017 % 0018 % Inputs : 0019 % y_in - input vector, to be projected 0020 % c - center of Lp ball 0021 % r - radius of Lp ball 0022 % p - exponent for Lp ball 0023 % 0024 % Outputs: 0025 % x_out - computed projection 0026 % n - number of iterations 0027 % 0028 % This is a matlab implementation, if possible the compiled c++ 0029 % implementation bpdq_proj_lpball_mex should be used instead. 0030 % This matlab file may be considered a more easily readable reference 0031 % implementation of the algorithm, and is included largely for 0032 % pedagogic purposes 0033 % 0034 % This file is part of BPDQ Toolbox (Basis Pursuit DeQuantizer) 0035 % Copyright (C) 2009, the BPDQ Team (see the file AUTHORS distributed with 0036 % this library) (See the notice at the end of the file.) 0037 0038 function [x_out,n] = bpdq_proj_lpball(y_in,c,r,p) 0039 tol=1e-15; 0040 maxiter=50; 0041 verbose=0; 0042 0043 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0044 % check first if input is inside specified ball 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 if norm(y_in - c,p)<r 0047 x_out=y_in; 0048 n=0; 0049 return 0050 end 0051 0052 if (p == Inf) 0053 x_out = c + sign(y_in - c).*min(abs(y_in - c), r); 0054 n=0; 0055 return 0056 end 0057 0058 if p<2 0059 error('Don''t try to use this code for p<2'); 0060 end 0061 0062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0063 % preprocessing to transform to c=0; r=1; case 0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0065 y=(y_in-c)/r; 0066 % preprocessing to remove signs 0067 % (work in orthant where everything is non-negative) 0068 ysign=sign(y); 0069 y=y.*ysign; 0070 0071 % now solve for xstar for c=0, r=1 0072 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0073 % Initialization 0074 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0075 N=numel(y); 0076 F=zeros(N+1,1); 0077 z=zeros(N+1,1); 0078 % xn1 : radial projection onto Lp ball 0079 xn1=y/norm(y,p); 0080 % xn2 : L\infty projection followed by radial projection onto Lp ball 0081 xn2=sign(y).*min(abs(y),1); 0082 xn2=xn2/norm(xn2,p); 0083 if (norm(xn1-y) < norm(xn2-y)) 0084 x0=xn1; 0085 else 0086 x0=xn2; 0087 end 0088 z(1:N)=x0; 0089 % pick least squares for lambda solving 0090 % lambda * x_i^(p-1) = y_i-x_i 0091 % as best it can 0092 lambda = x0.^(p-1) \ (y-x0); 0093 z(N+1)=lambda; 0094 0095 % Indices for constructing sparse Jacobian 0096 diag_i=1:N; 0097 diag_j=1:N; 0098 rightcol_i=1:N; 0099 rightcol_j=(N+1)*ones(1,N); 0100 bottomrow_i=(N+1)*ones(1,N); 0101 bottomrow_j=1:N; 0102 0103 all_i=[diag_i(:);rightcol_i(:);bottomrow_i(:)]; 0104 all_j=[diag_j(:);rightcol_j(:);bottomrow_j(:)]; 0105 0106 normz0=norm(z); % norm of input, useful for computing relative tolerance 0107 normF=tol+1; 0108 n=1; 0109 xn=z(1:N); 0110 % begin iteration 0111 0112 while normF/normz0 > tol; 0113 % Build residual F 0114 zpm1=z(1:N).^(p-1); % z to the p minus 1 0115 F(1:N)=z(1:N)+z(N+1)*zpm1-y; 0116 % divide by p to yield symmetric Jacobian ( ref my notes Nov 19) 0117 0118 F(N+1)= (sum(z(1:N).^p)-1)/p; 0119 normF=norm(F); 0120 0121 % Build Jacobian matrix J 0122 d=1+z(N+1)*(p-1)*z(1:N).^(p-2); % diagonal part 0123 0124 %% 0125 % This code computes the equivalent of inv(J)*-F using 0126 % the block inverse formula for J = [D,b ; b',0] 0127 % where D is diagonal. 0128 %% 0129 btwid=(zpm1./d); % inv(A)*b 0130 mu=zpm1'*btwid; % b^T*inv(A)*b 0131 chi=btwid'*(-F(1:N)); 0132 dz=[ [-F(1:N)./d+(-F(N+1)-chi)/mu*btwid]; (chi-(-F(N+1)) )/mu ]; 0133 0134 z=z+dz; 0135 0136 if verbose 0137 oxn=xn; 0138 xn=z(1:N); 0139 dxn=norm(xn - oxn)/norm(xn); 0140 cosn = dot(lpnormalvector(xn,p), (y- xn)/norm(y - xn)); 0141 fprintf('%i: |y-xn|_2=%e, |xn+1-xn|_2/|xn+1|=%e, cosn=%e, nF=%e\n',... 0142 n,norm(y - xn), dxn,cosn,normF); 0143 end 0144 n=n+1; 0145 if n>maxiter 0146 error('maximum # of iterations exceeded in proj_lpball_newton'); 0147 end 0148 end 0149 0150 xstar=z(1:N); 0151 % postprocessing to restore signs 0152 xstar=xstar.*ysign; 0153 0154 % postprocessing to transform back to c!=0; r!=1 case 0155 x_out=r*xstar+c; 0156 0157 0158 % The BPDQ Toolbox is free software: you can redistribute it and/or modify 0159 % it under the terms of the GNU General Public License as published by 0160 % the Free Software Foundation, either version 3 of the License, or 0161 % (at your option) any later version. 0162 % 0163 % The BPDQ Toolbox is distributed in the hope that it will be useful, 0164 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0165 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0166 % GNU General Public License for more details. 0167 % 0168 % You should have received a copy of the GNU General Public License 0169 % along with The BPDQ Toolbox. If not, see <http://www.gnu.org/licenses/>.