regtools: tikhonov(U,s,V,b,lambda,x_0) - File Exchange - MATLAB Central
function [fnew,cpred,wpred] = subsearch(f0,res,fdef,Afunc,Lfunc,lambda,S)
% Searches in subspace spanned by columns of S for the optimal solution
% of the regularization problem
% (lambda)^2*\vert\vert L*(f-fdef)\vert\vert^2 + \vert\vert d-A*f\vert\vert^2
% f0 = Initial guess at location of the minimum
% res = d - A*f0, the current residual
% fdef = The default image
% Afunc = Name of function which applies A to a vector
% Lfunc = Name of function which applies L to a vector
% lambda = Weighting between solution seminorm and data misfit
% S = matrix with search directions as its columns
%
% fnew = Position of minimum within subspace
% cpred = Predicted value of \vert\vert d-A*fnew\vert\vert^2
% wpred = Predicted value of \vert\vert L*(fnew-fdef)\vert\vert^2
% Sze Tan, University of Auckland, July 1998
nsrch = size(S,2);
pref = feval(Lfunc,f0-fdef);
w0 = pref’ * pref;
c0 = res’ * res;
fprintf(’Square of regularization semi-norm (current) = %f\n’,w0);
fprintf(’Square of data misfit norm (current) = %f\n’,c0);
AS = zeros(length(res),nsrch);
LS = zeros(length(pref),nsrch);
for k = 1 : nsrch
AS
,k) = feval(Afunc,S
,k));
LS
,k) = feval(Lfunc,S
,k));
end
Hc = AS’ * AS; Hw = LS’ * LS;
Gc = AS’ * res; Gw = -LS’*pref;
c = (Hc + lambda^2 * Hw) \ (Gc + lambda^2 *Gw);
cpred = c0 + c’*Hc*c - c’*Gc - Gc’*c;
wpred = w0 + c’*Hw*c - c’*Gw - Gw’*c;
fprintf(’Square of regularization semi-norm (predicted) = %f\n’,wpred);
fprintf(’Square of data misfit norm (predicted) = %f\n\n’,cpred);
fnew = f0 + S * c;