L-Curve算法拐点判据与最大曲率点搜索问题
本帖最后由 mxlzhenzhu 于 2013-11-13 23:20 编辑最小二乘中有一种L-Curve算法，此为背景。
现在得到一个离散点变化趋势，如下图1所示。
图1 离散点变化趋势
数据特征是，
1)，在该图左下角有一个曲率最大的点P，目前看起来有两个点(A&B)距离这个真实的曲率最大点也许很近了，但还不知道有多远，假如这个距离是distance。
2)，可以计算这条曲线的曲率，因为有曲率表达式；但是不可以对表达式求导。在【0,1】区间，曲率变化理论趋势为图2，图3是图1中所有点的曲率，包括A&B，这是因为这个distance还是太大了，漏掉了尖峰，就算有了尖峰还不一定就是曲率k取最大值的点：
图2 L曲线的曲率理论变化趋势，那个尖峰所在区间可能很小图3 如果分辨率不够就漏掉尖峰了
3），假如二分法去搜索这个点，会让计算量不可接受。
求算法。
本帖最后由 mxlzhenzhu 于 2013-11-13 23:45 编辑
他这个函数不全：http://cda.psych.uiuc.edu/html/4regutools/l_corner.html
还是知识产权问题啊。
看完了，算法有点复杂，自己慢慢搞吧。
http://www.ncbi.nlm.nih.gov/pubmed/11008433 过来学习了参考了 {:{39}:} 寂寞的部落 发表于 2013-11-28 16:00 static/image/common/back.gif
function = l_curve(U,sm,b,method,L,V)
%L_CURVE Plot the L-curve and find its "corner".
%
% =
% l_curve(U,s,b,method)
% l_curve(U,sm,b,method),sm =
% l_curve(U,s,b,method,L,V)
%
% Plots the L-shaped curve of eta, the solution norm || x || or
% semi-norm || L x ||, as a function of rho, the residual norm
% || A x - b ||, for the following methods:
% method = 'Tikh': Tikhonov regularization (solid line )
% method = 'tsvd': truncated SVD or GSVD (o markers)
% method = 'dsvd': damped SVD or GSVD (dotted line)
% method = 'mtsvd' : modified TSVD (x markers)
% The corresponding reg. parameters are returned in reg_param.If no
% method is specified then 'Tikh' is default.For other methods use plot_lc.
%
% Note that 'Tikh', 'tsvd' and 'dsvd' require either U and s (standard-
% form regularization) or U and sm (general-form regularization), while
% 'mtvsd' requires U and s as well as L and V.
%
% If any output arguments are specified, then the corner of the L-curve
% is identified and the corresponding reg. parameter reg_corner is
% returned.Use routine l_corner if an upper bound on eta is required.
% Reference: P. C. Hansen & D. P. O'Leary, "The use of the L-curve in
% the regularization of discrete ill-posed problems",SIAM J. Sci.
% Comput. 14 (1993), pp. 1487-1503.
% Per Christian Hansen, IMM, July 26, 2007.
% Set defaults.
if (nargin==3), method='Tikh'; end% Tikhonov reg. is default.
npoints = 200;% Number of points on the L-curve for Tikh and dsvd.
smin_ratio = 16*eps;% Smallest regularization parameter.
% Initialization.
= size(U); = size(sm);
if (nargout > 0), locate = 1; else locate = 0; end
beta = U'*b; beta2 = norm(b)^2 - norm(beta)^2;
if (ps==1)
s = sm; beta = beta(1:p);
else
s = sm(p:-1:1,1)./sm(p:-1:1,2); beta = beta(p:-1:1);
end
xi = beta(1:p)./s;
if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))
eta = zeros(npoints,1); rho = eta; reg_param = eta; s2 = s.^2;
reg_param(npoints) = max();
ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
for i=1:npoints
f = s2./(s2 + reg_param(i)^2);
eta(i) = norm(f.*xi);
rho(i) = norm((1-f).*beta(1:p));
end
if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
marker = '-'; txt = 'Tikh.';
elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4))
eta = zeros(p,1); rho = eta;
eta(1) = abs(xi(1))^2;
for k=2:p, eta(k) = eta(k-1) + abs(xi(k))^2; end
eta = sqrt(eta);
if (m > n)
if (beta2 > 0), rho(p) = beta2; else rho(p) = eps^2; end
else
rho(p) = eps^2;
end
for k=p-1:-1:1, rho(k) = rho(k+1) + abs(beta(k+1))^2; end
rho = sqrt(rho);
reg_param = (1:p)'; marker = 'o';
if (ps==1)
U = U(:,1:p); txt = 'TSVD';
else
U = U(:,1:p); txt = 'TGSVD';
end
elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))
eta = zeros(npoints,1); rho = eta; reg_param = eta;
reg_param(npoints) = max();
ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
for i=1:npoints
f = s./(s + reg_param(i));
eta(i) = norm(f.*xi);
rho(i) = norm((1-f).*beta(1:p));
end
if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
marker = ':';
if (ps==1), txt = 'DSVD'; else txt = 'DGSVD'; end
elseif (strncmp(method,'mtsv',4))
if (nargin~=6)
error('The matrices L and V must also be specified')
end
= size(L); rho = zeros(p,1); eta = rho;
= qr(L*V(:,n:-1:n-p),0);
for i=1:p
k = n-p+i;
Lxk = L*V(:,1:k)*xi(1:k);
zk = R(1:n-k,1:n-k)\(Q(:,1:n-k)'*Lxk); zk = zk(n-k:-1:1);
eta(i) = norm(Q(:,n-k+1:p)'*Lxk);
if (i < p)
rho(i) = norm(beta(k+1:n) + s(k+1:n).*zk);
else
rho(i) = eps;
end
end
if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
reg_param = (n-p+1:n)'; txt = 'MTSVD';
U = U(:,reg_param); sm = sm(reg_param);
marker = 'x'; ps = 2;% General form regularization.
else
error('Illegal method')
end
% Locate the "corner" of the L-curve, if required.
if (locate)
= l_corner(rho,eta,reg_param,U,sm,b,method);
end
% Make plot.
plot_lc(rho,eta,marker,ps,reg_param);
if locate
ax = axis;
HoldState = ishold; hold on;
loglog(,,':r',...
,,':r')
title(['L-curve, ',txt,' corner at ',num2str(reg_corner)]);
axis(ax)
if (~HoldState), hold off; end
end
function = mtsvd(U,s,V,b,k,L)
%MTSVD Modified truncated SVD regularization.
%
% = mtsvd(U,s,V,b,k,L)
%
% Computes the modified TSVD solution:
% x_k = V*[ xi_k ] .
% [ xi_0 ]
% Here, xi_k defines the usual TSVD solution
% xi_k = inv(diag(s(1:k)))*U(:,1:k)'*b ,
% and xi_0 is chosen so as to minimize the seminorm || L x_k ||.
% This leads to choosing xi_0 as follows:
% xi_0 = -pinv(L*V(:,k+1:n))*L*V(:,1:k)*xi_k .
%
% The truncation parameter must satisfy k > n-p.
%
% If k is a vector, then x_k is a matrix such that
% x_k = [ x_k(1), x_k(2), ... ] .
%
% The solution and residual norms are returned in eta and rho.
% Reference: P. C. Hansen, T. Sekii & H. Shibahashi, "The modified
% truncated-SVD method for regularization in general form", SIAM J.
% Sci. Stat. Comput. 13 (1992), 1142-1150.
% Per Christian Hansen, IMM, 12/22/95.
% Initialization.
m = size(U,1); = size(L);
lk = length(k); kmin = min(k);
if (kmin<n-p+1 | max(k)>n)
error('Illegal truncation parameter k')
end
x_k = zeros(n,lk);
beta = U(:,1:n)'*b; xi = beta./s;
eta = zeros(lk,1); rho =zeros(lk,1);
% Compute large enough QR factorization.
= qr(L*V(:,n:-1:kmin+1),0);
% Treat each k separately.
for j=1:lk
kj = k(j); xtsvd = V(:,1:kj)*xi(1:kj);
if (kj==n)
x_k(:,j) = xtsvd;
else
z = R(1:n-kj,1:n-kj)\(Q(:,1:n-kj)'*(L*xtsvd));
z = z(n-kj:-1:1);
x_k(:,j) = xtsvd - V(:,kj+1:n)*z;
end
eta(j) = norm(x_k(:,j));
rho(j) = norm(beta(kj+1:n) + s(kj+1:n).*z);
end
if (nargout > 1 & m > n)
rho = sqrt(rho.^2 + norm(b - U(:,1:n)*beta)^2);
end
function = tsvd(U,s,V,b,k)
%TSVD Truncated SVD regularization.
%
% = tsvd(U,s,V,b,k)
%
% Computes the truncated SVD solution
% x_k = V(:,1:k)*inv(diag(s(1:k)))*U(:,1:k)'*b .
% If k is a vector, then x_k is a matrix such that
% x_k = [ x_k(1), x_k(2), ... ] .
%
% The solution and residual norms are returned in eta and rho.
% Per Christian Hansen, IMM, 12/21/97.
% Initialization.
= size(V); lk = length(k);
if (min(k)<0 | max(k)>p)
error('Illegal truncation parameter k')
end
x_k = zeros(n,lk);
eta = zeros(lk,1); rho = zeros(lk,1);
beta = U(:,1:p)'*b;
xi = beta./s;
% Treat each k separately.
for j=1:lk
i = k(j);
if (i>0)
x_k(:,j) = V(:,1:i)*xi(1:i);
eta(j) = norm(xi(1:i));
rho(j) = norm(beta(i+1:p));
end
end
if (nargout > 1 & size(U,1) > p)
rho = sqrt(rho.^2 + norm(b - U(:,1:p)*beta)^2);
end
页:
[1]