FISTA_secondOrder.m

function [F, nr_of_iterations, Res] = FISTA_secondOrder(G, sampling, kappa, P_OmegaC, z, kmax)

%% FISTA_SECONDORDER: Solves the l1xl1 minimization problem for the splitting and/ or completion problem
% in Born approximation of order 2 numerically by using the FISTA algorithm.
%
% INPUT:    G           Observed far field matrix, nxhat*nd-array.
%           sampling    Structure containing information about the discretization.
%           kappa       Wave number, >0.
%           P_OmegaC    Projection operator corresponding to observable part Omega^C, nxhat*nd-array
%                       with only 0- and 1-entries, only 1-entries for splitting only.
%           z           Positions of the individual scatterers, 2*nr_of_sources-array,
%                       2*1-array for completion only.
%           kmax        Maximal number of FISTA iterations for stopping criterion.
%
% OUTPUT:   F                   Structure containing numerical reconstructions for solutions of
%                               l1xl1 minimization problem.
%           nr_of_iterations    Number of performed FISTA iterations.
%           Res                 Residuum-norms for all FISTA iterations, vector of length nr_of_iterations.
%
% SYNTAX: FISTA_secondOrder(G, sampling, kappa, P_OmegaC, z, kmax)
%
% *********************************************************************************************************

nr_of_sources = size(z,2);

% split and complete far field using iterated soft shrinkage:
landweberparam.kmax = kmax;  % max nr of iterations of iterated soft shrinkage algorithm
landweberparam.tol = 1e-3; % tolerance to be met in the Landweber iteration (stopping criterion)
landweberparam.omega = 1/(nr_of_sources+1); % relaxation parameter in the Landweber itertation
softshrinkparam.mu = 1e-3; % threshold parameter in the soft shrinkage operatore
softshrinkparam.type = 'complex'; % type of soft-shrinkage to be used ('real'/'complex')

softshrinkparam.weights = landweberparam.omega * softshrinkparam.mu;

% fast iteraitve thresholding algorithm:

[F, nr_of_iterations, Res] = fista_secondOrder(G, sampling, kappa, P_OmegaC, z, landweberparam, softshrinkparam);

indOfInterest = sub2ind([nr_of_sources,nr_of_sources],1:nr_of_sources,1:nr_of_sources);

F =  {mysum(F, nr_of_sources^2), F{indOfInterest}};

end


function [F, nr_of_iterations, Res] = fista_secondOrder(G, sampling, kappa, P_OmegaC, z, landweberparam, softshrinkparam) 

% initialize some variables:
nxhat = sampling.nxhat;
nd = sampling.nd;

nr_of_sources = size(z,2);

[X,Y] = meshgrid(1:nr_of_sources, 1:nr_of_sources);
indices = [reshape(X,nr_of_sources^2,1),reshape(Y,nr_of_sources^2,1)]; % (J^2+1)*2-array that indicates the finite dimensional subspace where the corresponding block should lie in.
% (0,0) indicates $V_\Omega$; (j,l) indicates $V_{N_l,N_j}$, i.e.
% cutoff-parameter N_l in row direction und N_j in column direction
clear X Y

% initialize some more variables:
X = cell(nr_of_sources^2,1);  % columns of B correspond to the vectors b_j in the manuscript
for iters=1:nr_of_sources^2
    X{iters} = zeros(nxhat,nd);
end
KstarRes = cell(nr_of_sources^2,1);
iteri = 1;  % iteration counter
Res = zeros(landweberparam.kmax, 1);  % norm of the residual in each step of the iteration
Res(1) = landweberparam.tol+1;  % just needed to get iteration started; will be removed later

Y = X; Xalt = X;
Yhelp = Y;
talt = 1;

% disp('FISTA Algorithm starts.')
% tic

% fast iterated soft shrinkage:
while (Res(iteri) > landweberparam.tol) && (iteri < landweberparam.kmax)
   
    % apply operator K:
    for iters = 1 : nr_of_sources^2
       incidentDir = indices(iters,1);
       detectorPos = indices(iters,2);
       Yhelp{iters} = fft2(translOp(ifft2(Y{iters})*(nxhat*nd),sampling,-z(:,incidentDir),-z(:,detectorPos), kappa))/(nxhat*nd);
       clear incidentDir detectorPos
    end    
    KB = ifft2(mysum(Yhelp,nr_of_sources^2))*(nxhat*nd);
    KB = KB .* P_OmegaC;  % apply P_OmegaC
    
    % evaluate residual:
    Res = G - KB;
    
    % apply operator K^*:
    Res = Res .* P_OmegaC;  % apply P_OmegaC
    Res = fft2(Res)/(nxhat*nd);
    for iters = 1 : nr_of_sources^2
       incidentDir = indices(iters,1);
       detectorPos = indices(iters,2);
       KstarRes{iters} = fft2(translOp(ifft2(Res)*(nxhat*nd),sampling,z(:,incidentDir),z(:,detectorPos),kappa))/(nxhat*nd);
       clear incidentDir detectorPos
    end
    
    X = myplus(Y, myprod(landweberparam.omega,KstarRes,nr_of_sources^2),nr_of_sources^2);
    
    % apply soft shrinkage:
    X = softshrink(X, softshrinkparam.weights, softshrinkparam.type, nr_of_sources^2);
    
    t = (1+sqrt(1+4*talt^2))/2;
    
    Y = myplus(X,myprod((talt-1)/t,(myplus(X,myprod(-1,Xalt,nr_of_sources^2),nr_of_sources^2)),nr_of_sources^2),nr_of_sources^2);
    
    Xalt = X;
    talt = t;
    
    Res(iteri+1) = norm(Res);  % HS norm of residual    
    iteri = iteri+1;  % update iteration counter
    
end

% toc
% disp('FISTA Algorithm finished.')
% 
% fprintf('\n\nNumber of iterations: %i\n', iteri);

nr_of_iterations = iteri;

% shift back and compute inverse Fourier transform:

F = cell(nr_of_sources^2, 1);  % cell array containing splitted far field components
for iters = 1 : nr_of_sources^2
    Fhelp = ifft2(X{iters})*(nxhat*nd);
    incidentDir = indices(iters,1);
    detectorPos = indices(iters,2);
    F{iters} = translOp(Fhelp,sampling,-z(:,incidentDir),-z(:,detectorPos),kappa);
    clear incidentDir detectorPos
end


end


function Bsum = mysum(B, nr_of_sources)

Bsum = B{1};
for iterk = 2:nr_of_sources
    
    Bsum = Bsum + B{iterk};
    
end
end


function Z = myplus(X, Y, nr_of_sources)

Z = cell(nr_of_sources, 1);
for iterk = 1:nr_of_sources
    
    Z{iterk} = X{iterk} + Y{iterk};
    
end
end


function Z = myprod(alpha, X, nr_of_sources)

Z = cell(nr_of_sources, 1);
for iterk = 1:nr_of_sources
    
    Z{iterk} = alpha*X{iterk};
    
end
end