Update FISTA_secondOrder.m

src/FISTA_secondOrder.m

-function [F, nr_of_iterations, Resnorm] = FISTA_secondOrder(G, sampling, kappa, P_OmegaC, z, kmax)
-%% L1SCHEME:
-%  INPUT:   farfield - far field matrix
-%           kappa - wave number
-%           Omega - part of [0,2pi] where the far field pattern is unkown
-%           z - 2xm-vector containing the centers of the circles describing 
-%               the a priori information on the source locations
-%           r - 1xm-vector containing the radii of the circles describing 
-%               the a priori information on the source locations
-%  OUTPUT:  f - cell array containing the m source components
-%  EXAMPLE USAGE: 
-
-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, Resnorm] = 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, Resnorm] = fista_secondOrder(G, sampling, kappa, P_OmegaC, z, landweberparam, softshrinkparam)
-%% FISTA:
-%  INPUT:   umeas - far field pattern
-%           kappa - wave number
-%           Omega - 
-%           z - 2xnr_of_sources vector containing presumed source positions
-%           landweberparam.Nmax -- max nr of iterations
-%           landweberparam.tol -- tolerance to be met in the Landweber iteration
-%           landweberparam.omega -- relaxation parameter in the Landweber itertation
-%           softshrinkparam.mu -- threshold parameter in the soft shrinkage operatore
-%           softshrinkparam.type -- type of soft-shrinkage to be used ('real'/'complex')
-%  OUTPUT:  f - cell array containing split far field components
-%           
-%  EXAMPLE USAGE: 
-
-
-% 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
-Resnorm = zeros(landweberparam.kmax, 1);  % norm of the residual in each step of the iteration
-Resnorm(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 (Resnorm(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;
-    
-    Resnorm(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
\ No newline at end of file
+function [F, nr_of_iterations, Res] = FISTA_secondOrder(G, sampling, kappa, P_OmegaC, z, kmax)
+
+% 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.
+%
+% *********************************************************************************************************
+
+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