function [gd,relres,iter] = gabconvexopt(g,a,M,varargin)
%GABCONVEXOPT Compute a window using convex optimization
% Usage: gout=gabconvexopt(g,a,M);
% gout=gabconvexopt(g,a,M, varagin);
%
% Input parameters:
% g : Window function /initial point (tight case)
% a : Time shift
% M : Number of Channels
%
% Output parameters:
% gout : Output window
% iter : Number of iterations
% relres : Reconstruction error
%
% GABCONVEXOPT(g,a,M) computes a window gout which is the optimal
% solution of the convex optimization problem below
%
% gd = argmin_x || alpha x||_1 + || beta Fx||_1
%
% + || omega (x -g_l) ||_2^2 + delta || x ||_S0
%
% + gamma || nabla F x ||_2^2 + mu || nabla x ||_2^2
%
% such that x satifies the constraints
%
% Three constraints are possible:
%
% x is dual with respect of g
%
% x is tight
%
% x is compactly supported on Ldual
%
% *Note**: This function require the unlocbox. You can download it at
% http://unlocbox.sourceforge.net
%
% The function uses an iterative algorithm to compute the approximate.
% The algorithm can be controlled by the following flags:
%
% 'alpha',alpha Weight in time. If it is a scalar, it represent the
% weights of the entire L1 function in time. If it is a
% vector, it is the associated weight assotiated to each
% component of the L1 norm (length: Ldual).
% Default value is alpha=0.
% *Warning**: this value should not be too big in order to
% avoid the the L1 norm proximal operator kill the signal.
% No L1-time constraint: alpha=0
%
% 'beta',beta Weight in frequency. If it is a scalar, it represent the
% weights of the entire L1 function in frequency. If it is a
% vector, it is the associated weight assotiated to each
% component of the L1 norm in frequency. (length: Ldual).
% Default value is beta=0.
% *Warning**: this value should not be too big in order to
% avoid the the L1 norm proximal operator kill the signal.
% No L1-frequency constraint: beta=0
%
% 'omega',omega Weight in time of the L2-norm. If it is a scalar, it represent the
% weights of the entire L2 function in time. If it is a
% vector, it is the associated weight assotiated to each
% component of the L2 norm (length: Ldual).
% Default value is omega=0.
% No L2-time constraint: omega=0
%
% 'glike',g_l g_l is a windows in time. The algorithm try to shape
% the dual window like g_l. Normalization of g_l is done
% automatically. To use option omega should be different
% from 0. By default g_d=0.
%
% 'mu', mu Weight of the smooth constraint Default value is 1.
% No smooth constraint: mu=0
%
% 'gamma', gamma Weight of the smooth constraint in frequency. Default value is 1.
% No smooth constraint: gamma=0
%
% 'delta', delta Weight of the S0-norm. Default value is 0.
% No S0-norm: delta=0
%
% 'support' Ldual Add a constraint on the support. The windows should
% be compactly supported on Ldual.
%
% 'tight' Look for a tight windows
%
% 'dual' Look for a dual windows (default)
%
% 'painless' Construct a starting guess using a painless-case
% approximation. This is the default
%
% 'zero' Choose a starting guess of zero.
%
% 'rand' Choose a random starting phase.
%
% 'tol',t Stop if relative residual error is less than the
% specified tolerance.
%
% 'maxit',n Do at most n iterations. default 200
%
% 'print' Display the progress.
%
% 'debug' Display all the progresses.
%
% 'quiet' Don't print anything, this is the default.
%
% 'fast' Fast algorithm, this is the default.
%
% 'slow' Safer algorithm, you can try this if the fast algorithm
% is not working. Before using this, try to iterate more.
%
% 'printstep',p If 'print' is specified, then print every p'th
% iteration. Default value is p=10;
%
% 'hardconstraint' Force the projection at the end (default)
%
% 'softconstaint' Do not force the projection at the end
%
% See also: gaboptdual, gabdual, gabtight, gabfirtight, gabopttight
%
% Url: http://ltfat.github.io/doc/gabor/gabconvexopt.html
% Copyright (C) 2005-2023 Peter L. Soendergaard <peter@sonderport.dk> and others.
% This file is part of LTFAT version 2.6.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% Author: Nathanael Perraudin
% Date : 18 Feb 2014
if nargin<4
error('%s: Too few input parameters.',upper(mfilename));
end;
if numel(g)==1
error('g must be a vector (you probably forgot to supply the window function as input parameter.)');
end;
definput.keyvals.L=[];
definput.keyvals.lt=[0 1];
definput.keyvals.tol=1e-6;
definput.keyvals.maxit=200;
definput.keyvals.printstep=10;
definput.flags.print={'quiet','print','debug'};
definput.flags.algo={'fast','slow'};
definput.flags.constraint={'hardconstraint','softconstaint'};
definput.flags.startphase={'painless','zero','rand'};
definput.flags.type={'dual','tight'};
definput.keyvals.alpha=0;
definput.keyvals.omega=0;
definput.keyvals.beta=0;
definput.keyvals.mu=1;
definput.keyvals.gamma=1;
definput.keyvals.vart=0;
definput.keyvals.varf=0;
definput.keyvals.var2t=0;
definput.keyvals.var2f=0;
definput.keyvals.support=0;
definput.keyvals.delta=0;
definput.keyvals.deltaw=0;
definput.keyvals.glike=zeros(size(g));
[flags,kv]=ltfatarghelper({'L','tol','maxit'},definput,varargin);
% Determine the window. The window /must/ be an FIR window, so it is
% perfectly legal to specify L=[] when calling gabwin
[g,info]=gabwin(g,a,M,[],kv.lt,'callfun',upper(mfilename));
if kv.support
Ldual=kv.support;
% Determine L. L must be longer than L+Ldual+1 to make sure that no convolutions are periodic
L=dgtlength(info.gl+Ldual+1,a,M);
else
L=length(g);
Ldual=L;
end
b=L/M;
% Determine the initial guess
if flags.do_zero
gd_initial=zeros(Ldual,1);
end;
if flags.do_rand
gd_initial=rand(size(g));
end;
if flags.do_painless
gsmall=long2fir(g,M);
gdsmall=gabdual(gsmall,a,M);
gd_initial=fir2long(gdsmall,Ldual);
end;
% -------- do the convex optimization stuff
% Define the long original window
glong=fir2long(g,L);
%gabframebounds(g,a,M)
% Initial point
xin=gd_initial;
xin=fir2long(xin,L);
% -- * Setting the different prox for ppxa *--
% ppxa will minimize all different proxes
% value test for the selection constraint
nb_priors=0;
% - variance -
if kv.vart % constraint in time
if flags.do_debug
param_l1.verbose=1; % display the results
else
param_l1.verbose=0; % do not display anything
end
% alpha is a scalar
if mod(L,2)
w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
else
w=[0:1:L/2-1,L/2:-1:1]';
end
w=w.^2/L;
param_l1.weights=w;
nb_priors=nb_priors+1;
g11.prox= @(x,T) prox_l1(x,kv.vart*T,param_l1); % define the prox_l1 as operator
g11.eval= @(x) kv.vart*norm(w.*x,1); % the objectiv function is the l1 norm
else % no L1 in time constraint
g11.prox= @(x,T) x;
g11.eval= @(x) 0;
end
% - variance -
if kv.varf % constraint in time
param_l1_fourier.A= @(x) 1/sqrt(L)*fft(x); % Fourier operator
param_l1_fourier.At= @(x) sqrt(L)*ifft(x); % adjoint of the Fourier operator
if flags.do_debug
param_l1_fourier.verbose=1; % display the results
else
param_l1_fourier.verbose=0; % do not display anything
end
if mod(L,2)
w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
else
w=[0:1:L/2-1,L/2:-1:1]';
end
w=w.^2/L;
param_l1_fourier.weights=w;
nb_priors=nb_priors+1;
g12.prox= @(x,T) prox_l1(x,kv.varf*T,param_l1_fourier); % define the prox_l1 as operator
g12.eval= @(x) kv.varf*norm(w.*x,1); % the objectiv function is the l1 norm
else % no L1 in time constraint
g12.prox= @(x,T) x;
g12.eval= @(x) 0;
end
% - variance2 -
if kv.var2t % constraint in time
if flags.do_debug
param_l2.verbose=1; % display the results
else
param_l2.verbose=0; % do not display anything
end
% alpha is a scalar
if mod(L,2)
w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
else
w=[0:1:L/2-1,L/2:-1:1]';
end
w=w/sqrt(L);
param_l2.weights=w;
nb_priors=nb_priors+1;
g13.prox= @(x,T) prox_l2(x,kv.var2t*T,param_l2); % define the prox_l1 as operator
g13.eval= @(x) kv.var2t*norm(w.*x,2)^2; % the objectiv function is the l1 norm
else % no L1 in time constraint
g13.prox= @(x,T) x;
g13.eval= @(x) 0;
end
% - variance2 -
if kv.var2f % constraint in time
param_l2_fourier.A= @(x) 1/sqrt(L)*fft(x); % Fourier operator
param_l2_fourier.At= @(x) sqrt(L)*ifft(x); % adjoint of the Fourier operator
if flags.do_debug
param_l2_fourier.verbose=1; % display the results
else
param_l2_fourier.verbose=0; % do not display anything
end
if mod(L,2)
w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
else
w=[0:1:L/2-1,L/2:-1:1]';
end
w=w/sqrt(L);
param_l2_fourier.weights=w;
nb_priors=nb_priors+1;
g14.prox= @(x,T) prox_l2(x,kv.var2f*T,param_l2_fourier); % define the prox_l1 as operator
g14.eval= @(x) kv.var2f*norm(w.*x,2)^2; % the objectiv function is the l1 norm
else % no L1 in time constraint
g14.prox= @(x,T) x;
g14.eval= @(x) 0;
end
% - small L1 norm in coefficient domain -
if kv.alpha % constraint in time
if flags.do_debug
param_l1.verbose=1; % display the results
else
param_l1.verbose=0; % do not display anything
end
if length(kv.alpha)==1 % alpha is a scalar
kv.alpha=ones(size(xin))*kv.alpha;
end
param_l1.weights=kv.alpha;
nb_priors=nb_priors+1;
g1.prox= @(x,T) prox_l1(x,T,param_l1); % define the prox_l1 as operator
g1.eval= @(x) norm(kv.alpha.*x,1); % the objectiv function is the l1 norm
else % no L1 in time constraint
g1.prox= @(x,T) x;
g1.eval= @(x) 0;
end
% - small L1 norm in Fourier domain -
if kv.beta %frequency constraint
param_l1_fourier.A= @(x) 1/sqrt(L)*fft(x); % Fourier operator
param_l1_fourier.At= @(x) sqrt(L)*ifft(x); % adjoint of the Fourier operator
if flags.do_debug
param_l1_fourier.verbose=1; % display the results
else
param_l1_fourier.verbose=0; % Do not display anything
end
if length(kv.beta)==1 % alpha is a scalar
kv.beta=ones(size(xin))*kv.beta;
end
param_l1_fourier.weights=kv.beta;
% Here are the step for the prox
% 2) go into the Fourier domain (prox_l1)
% 3) soft thresholding (prox_l1)
% 4) back in the time domain (prox_l1)
nb_priors=nb_priors+1;
g3.prox= @(x,T) prox_l1(x,T,param_l1_fourier);
g3.eval= @(x) norm(kv.beta.*fft(x),1); % objectiv function
else % no L1 in frequency constraint
g3.prox= @(x,T) x;
g3.eval= @(x) 0; % objectiv function
end
% - DUAL OR TIGHT?-
if flags.do_tight
% tight windows
g2.prox= @(x,T) gabtight(x,a,M); % set the prox
g2.eval= @(x) norm(x-gabdual(x,a,M,L)); % objectiv function
else
% - projection on a B2 ball -
% Frame-type matrix of the adjoint lattice
%G=tfmat('dgt',glong,M,a);
Fal=frame('dgt',glong,M,a);
G=framematrix(Fal,L);
d=[a/M;zeros(a*b-1,1)];
% Using a B2 ball projection
% || Gcut' x - b ||_2 < epsilon
% param_proj.A = @(x) G'*x; % forward operator
% param_proj.At = @(x) G*x; % adjoint operator
% param_proj.y = d;
% param_proj.maxit = 200; % maximum of iteration
% param_proj.tight=0; % not a tight frame
% param_proj.nu=norm(G)^2; % frame bound on Gcut'
% param_proj.verbose=0; % diplay summary at the end
% param_proj.epsilon=10*eps; % radius of the B2 ball
% g2.prox= @(x,T) fast_proj_b2(x,T,param_proj); % set the prox
% Using a direct projection (better solution)
param_proj.verbose=flags.do_debug;
param_proj.y=d;
param_proj.A=G';
param_proj.AAtinv=(G'*G)^(-1);
g2.prox= @(x,T) proj_dual(x,T,param_proj); % set the prox
g2.eval= @(x) norm(G'*x-d); % objectiv function
end
% SUPPORT CONSTRAINT
if kv.support
% - set null coefficient
g4.prox = @(x,T) forceeven(fir2long(long2fir(x,Ldual),L));
g4.eval = @(x) 0;
% - function apply the two projections thanks to a poc algorithm.
if flags.do_tight
G={g2,g4};
paramPOCS.tol=20*eps;
paramPOCS.maxit=5000;
paramPOCS.verbose=flags.do_print+flags.do_debug;
paramPOCS.abs_tol=1;
g5.prox = @(x,T) pocs(x,G,paramPOCS);
% g5.prox = @(x,T) ppxa(x,G,paramPOCS);
% g5.prox = @(x,T) douglas_rachford(x,g2,g4,paramPOCS);
% g5.prox = @(x,T) pocs2(x,g2,g4,20*eps,2000, flags.do_print+flags.do_debug);
g5.eval = @(x) 0;
else
Fal=frame('dgt',glong,M,a);
G=framematrix(Fal,L);
d=[a/M;zeros(a*b-1,1)];
Lfirst=ceil(Ldual/2);
Llast=Ldual-Lfirst;
Gcut=G([1:Lfirst,L-Llast+1:L],:);
param_proj2.verbose=flags.do_debug;
param_proj2.y=d;
param_proj2.A=Gcut';
param_proj2.AAtinv=pinv(Gcut'*Gcut);
g5.prox= @(x,T) fir2long(proj_dual(long2fir(x,Ldual),T,param_proj2),L); % set the prox
g5.eval= @(x) norm(G'*x-d); % objectiv function
end
else
g4.prox= @(x,T) x;
g4.eval= @(x) 0; % objectiv function
g5=g2;
end
% - function apply the two projections thanks to a douglas rachford algorithm.
% param_douglas.verbose=1;
% param_douglas.abs_tol=1;
% param_douglas.maxit=2000;
% param_douglas.tol=20*eps;
% g6.prox = @(x,T) douglas_rachford(x,g2,g4,param_douglas);
% g6.eval = @(x) 0;
% - small gradient norm -
% this is the smoothing parameter
if kv.mu
if flags.do_debug
param_l2grad.verbose=1; % display the results
else
param_l2grad.verbose=0; % Do not display anything
end
nb_priors=nb_priors+1;
g7.prox = @(x,T) prox_l2grad(fir2long(x,L),kv.mu*T,param_l2grad);
g7.eval = @(x) norm(gradient(x))^2;
else
g7.prox = @(x,T) x;
g7.eval = @(x) 0;
end
% - small gradient norm in fourrier-
% this is the smoothing parameter
if kv.gamma
if flags.do_debug
param_l2grad.verbose=1; % display the results
else
param_l2grad.verbose=0; % Do not display anything
end
nb_priors=nb_priors+1;
g9.prox = @(x,T) prox_l2gradfourier(fir2long(x,L),kv.gamma*T,param_l2grad);
g9.eval = @(x) norm(gradient(1/sqrt(L)*fft(x)))^2;
else
g9.prox = @(x,T) x;
g9.eval = @(x) 0;
end
% - small L2 norm in coefficient domain -
if kv.omega % constraint in time
if flags.do_debug
param_l2.verbose=1; % display the results
else
param_l2.verbose=0; % do not display anything
end
if length(kv.omega)==1 % alpha is a scalar
kv.alpha=ones(size(xin))*kv.omega;
end
param_l2.weights=kv.omega;
if sum(kv.glike)
kv.glike=fir2long(kv.glike,L);
glike=kv.glike/norm(kv.glike)*norm(gabdual(g,a,M));
param_l2.y=fir2long(glike,L);
end
nb_priors=nb_priors+1;
g8.prox= @(x,T) prox_l2(x,T,param_l2); % define the prox_l2 as operator
g8.eval= @(x) norm(kv.omega.*x-kv.glike,'fro')^2; % the objectiv function is the l2 norm
else % no L1 in time constraint
g8.prox= @(x,T) x;
g8.eval= @(x) 0;
end
% - small S0 norm -
if kv.delta %frequency constraint
gauss=pgauss(L,1);
[A,B]=gabframebounds(gauss,1,L);
AB=(A+B)/2;
param_S0.A= @(x) dgt(x,gauss,1,L)/sqrt(AB);
param_S0.At= @(x) idgt(x,gauss,1,L)/sqrt(AB);
if flags.do_debug
param_S0.verbose=1; % display the results
else
param_S0.verbose=0; % Do not display anything
end
nb_priors=nb_priors+1;
g10.prox= @(x,T) prox_l1(x,T*kv.delta,param_S0);
g10.eval= @(x) kv.delta*norm(reshape(dgt(x,gauss,1,L),[],1),1); % objectiv function
else % no L1 in frequency constraint
g10.prox= @(x,T) x;
g10.eval= @(x) 0; % objectiv function
end
% - small weighted S0 norm -
if kv.deltaw %frequency constraint
gauss=pgauss(L,1);
[A,B]=gabframebounds(gauss,1,L);
AB=(A+B)/2;
param_S0.A= @(x) dgt(x,gauss,1,L)/sqrt(AB);
param_S0.At= @(x) idgt(x,gauss,1,L)/sqrt(AB);
if flags.do_debug
param_S0.verbose=1; % display the results
else
param_S0.verbose=0; % Do not display anything
end
if mod(L,2)
w=[0:1:(L-1)/2,(L-1)/2:-1:1]';
else
w=[0:1:L/2-1,L/2:-1:1]';
end
w=w/sqrt(L);
%W=w*w';
W=repmat(w,1,L).^2+repmat(w',L,1).^2;
W=sqrt(W);
param_S0.weights=W;
nb_priors=nb_priors+1;
g15.prox= @(x,T) prox_l1(x,T*kv.deltaw,param_S0);
g15.eval= @(x) kv.deltaw*norm(reshape(dgt(x,gauss,1,L),[],1),1); % objectiv function
else % no L1 in frequency constraint
g15.prox= @(x,T) x;
g15.eval= @(x) 0; % objectiv function
end
% -- * PPXA function, the solver * --
% parameter for the solver
param.maxit=kv.maxit; % maximum number of iteration
param.tol=kv.tol;
if flags.do_quiet
param.verbose=0;
end
% Definition of the function f (the order is important)
if flags.do_fast && flags.do_tight
F={g1, g3,g7,g9,g8, g2, g4,g10,g11,g12,g13,g14,g15};
else
F={g1, g3,g7,g9,g8, g5,g10,g11,g12,g13,g14,g15};
end
% solving the problem
if nb_priors
[gd,iter,~]=ppxa(xin,F,param);
% Force the hard constraint
if flags.do_hardconstraint
% In case of use of the douglas rachford algo instead of POCS
% gd=g6.prox(gd,0); % force the constraint
gd=g5.prox(gd,0);
end
else
fprintf( ' Warning!!! No prior selected! -- Only perform a projection. \n')
gd=g5.prox(xin,0);
end
% compute the error
if flags.do_tight
relres=gabdualnorm(gd,gd,a,M,L);
else
relres=gabdualnorm(g,gd,a,M,L);
end
if kv.support
% set the good size
gd=long2fir(gd,Ldual);
end
end
% function x=pocs2(x,g1,g2,tol,maxii,flagp)
% % this function implement a POCS algorithm, projection onto convex Set
% % using the differents projection of the algorithm.
% tola=1;
% ii=0;
% tola_old=tola;
% while (tola>tol)
% x=g2.prox(g1.prox(x,0),0);
% tola=g1.eval(x);
% ii=ii+1;
% if (logical(1-logical(mod(ii,50))) && flagp)
% fprintf(' POCS sub-iteration: %i -- Tol : %g\n',ii,tola)
% end
% if ii> maxii
% break;
% end
% if abs(tola_old-tola)/tola<tol % avoid infinite loop
% break;
% end
% tola_old=tola;
% end
% end
function x=forceeven(x)
% this function force the signal to be even
x= (x+involute(x))/2;
end