tikhonov.zip_L曲线正则化_Tikhonov 正则化_Tikhonov正则化_tikhonov_正则化

function [reg_c,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,s,b,method,M) %L_CORNER Locate the "corner" of the L-curve. % % [reg_c,rho_c,eta_c] = % l_corner(rho,eta,reg_param) % l_corner(rho,eta,reg_param,U,s,b,method,M) % l_corner(rho,eta,reg_param,U,sm,b,method,M) , sm = [sigma,mu] % % Locates the "corner" of the L-curve in log-log scale. % % It is assumed that corresponding values of || A x - b ||, || L x ||, % and the regularization parameter are stored in the arrays rho, eta, % and reg_param, respectively (such as the output from routine l_curve). % % If nargin = 3, then no particular method is assumed, and if % nargin = 2 then it is issumed that reg_param = 1:length(rho). % % If nargin >= 6, then the following methods are allowed: % method = 'Tikh' : Tikhonov regularization % method = 'tsvd' : truncated SVD or GSVD % method = 'dsvd' : damped SVD or GSVD % method = 'mtsvd' : modified TSVD, % and if no method is specified, 'Tikh' is default. If the Spline Toolbox % is not available, then only 'Tikh' and 'dsvd' can be used. % % An eighth argument M specifies an upper bound for eta, below which % the corner should be found. % The following functions from the Spline Toolbox are needed if % method differs from 'Tikh' or 'dsvd': % fnder, ppbrk, ppmak, ppual, sp2pp, sorted, spbrk, spmak, sprpp. % Per Christian Hansen, IMM, Dec. 12, 2002. % Set default regularization method. if (nargin <= 3) method = 'none'; if (nargin==2), reg_param = [1:length(rho)]'; end else if (nargin==6), method = 'Tikh'; end end % Set threshold for skipping very small singular values in the % analysis of a discrete L-curve. s_thr = eps; % Neglect singular values less than s_thr. % Set default parameters for treatment of discrete L-curve. deg = 2; % Degree of local smooting polynomial. q = 2; % Half-width of local smoothing interval. orig = 2 order = 4; % Order of fitting 2-D spline curve. % Initialization. if (length(rho) < order) error('Too few data points for L-curve analysis') end if (nargin > 3) [p,ps] = size(s); [m,n] = size(U); beta = U'*b; if (m>n), b0 = b - U*beta; end if (ps==2) s = s(p:-1:1,1)./s(p:-1:1,2); U = U(:,p:-1:1); beta = beta(p:-1:1); end xi = beta./s; end % Restrict the analysis of the L-curve according to M (if specified). if (nargin==8) index = find(eta < M); rho = rho(index); eta = eta(index); reg_param = reg_param(index); end if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4)) % The L-curve is differentiable; computation of curvature in % log-log scale is easy. % Compute g = - curvature of L-curve. g = lcfun(reg_param,s,beta,xi); % Locate the corner. If the curvature is negative everywhere, % then define the leftmost point of the L-curve as the corner. [gmin,gi] = min(g); reg_c = fminbnd('lcfun',... reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),... optimset('Display','off'),s,beta,xi); % Minimizer. kappa_max = - lcfun(reg_c,s,beta,xi); % Maximum curvature. if (kappa_max < 0) lr = length(rho); reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr); else f = (s.^2)./(s.^2 + reg_c^2); eta_c = norm(f.*xi); rho_c = norm((1-f).*beta); if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end end elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4) | ... strncmp(method,'mtsv',4) | strncmp(method,'none',4)) % The L-curve is discrete and may include unwanted fine-grained % corners. Use local smoothing, followed by fitting a 2-D spline % curve to the smoothed discrete L-curve. % Check if the Spline Toolbox exists, otherwise return. if (exist('splines')~=7) error('The Spline Toolbox in not available so l_corner cannot be used') end % For TSVD, TGSVD, and MTSVD, restrict the analysis of the L-curve % according to s_thr. if (nargin > 3) if (nargin==8) % In case the bound M is in action. s = s(index,:); beta = beta(index); xi = xi(index); end index = find(s > s_thr); rho = rho(index); eta = eta(index); reg_param = reg_param(index); s = s(index); beta = beta(index); xi = xi(index); end % Convert to logarithms. lr = length(rho); lrho = log(rho); leta = log(eta); slrho = lrho; sleta = leta; % For all interior points k = q+1:length(rho)-q-1 on the discrete % L-curve, perform local smoothing with a polynomial of degree deg % to the points k-q:k+q. v = [-q:q]'; A = zeros(2*q+1,deg+1); A(:,1) = ones(length(v),1); for j = 2:deg+1, A(:,j) = A(:,j-1).*v; end for k = q+1:lr-q-1 cr = A\lrho(k+v); slrho(k) = cr(1); ce = A\leta(k+v); sleta(k) = ce(1); end % Fit a 2-D spline curve to the smoothed discrete L-curve. sp = spmak([1:lr+order],[slrho';sleta']); pp = ppbrk(sp2pp(sp),[4,lr+1]); % Extract abscissa and ordinate splines and differentiate them. % Compute as many function values as default in spleval. P = spleval(pp); dpp = fnder(pp); D = spleval(dpp); ddpp = fnder(pp,2); DD = spleval(ddpp); ppx = P(1,:); ppy = P(2,:); dppx = D(1,:); dppy = D(2,:); ddppx = DD(1,:); ddppy = DD(2,:); % Compute the corner of the discretized .spline curve via max. curvature. % No need to refine this corner, since the final regularization % parameter is discrete anyway. % Define curvature = 0 where both dppx and dppy are zero. k1 = dppx.*ddppy - ddppx.*dppy; k2 = (dppx.^2 + dppy.^2).^(1.5); I_nz = find(k2 ~= 0); kappa = zeros(1,length(dppx)); kappa(I_nz) = -k1(I_nz)./k2(I_nz); [kmax,ikmax] = max(kappa); x_corner = ppx(ikmax); y_corner = ppy(ikmax); % Locate the point on the discrete L-curve which is closest to the % corner of the spline curve. Prefer a point below and to the % left of the corner. If the curvature is negative everywhere, % then define the leftmost point of the L-curve as the corner. if (kmax < 0) reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr); else index = find(lrho < x_corner & leta < y_corner); if (length(index) > 0) [dummy,rpi] = min((lrho(index)-x_corner).^2 + (leta(index)-y_corner).^2); rpi = index(rpi); else [dummy,rpi] = min((lrho-x_corner).^2 + (leta-y_corner).^2); end reg_c = reg_param(rpi); rho_c = rho(rpi); eta_c = eta(rpi); end elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4)) % The L-curve is differentiable; computation of curvature in % log-log scale is easy. % Compute g = - curvature of L-curve. g = lcfun(reg_param,s,beta,xi,1); % Locate the corner. If the curvature is negative everywhere, % then define the leftmost point of the L-curve as the corner. [gmin,gi] = min(g); reg_c = fminbnd('lcfun',... reg_param(min(gi+1,length(g))),reg_param(max(gi-1,1)),... optimset('Display','off'),s,beta,xi,1); % Minimizer. kappa_max = - lcfun(reg_c,s,beta,xi,1); % Maximum curvature. if (kappa_max < 0) lr = length(rho); reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr); else f = s./(s + reg_c); eta_c = norm(f.*xi); rho_c = norm((1-f).*beta); if (m>n), rho_c = sqrt(rho_c^2 + norm(b0)^2); end end else, error('Illegal method'), end