L-Curve算法拐点判据与最大曲率点搜索问题

mxlzhenzhu

最小二乘中有一种L-Curve算法，此为背景。

现在得到一个离散点变化趋势，如下图1所示。

图1 离散点变化趋势

数据特征是，
1)，在该图左下角有一个曲率最大的点P，目前看起来有两个点(A&B)距离这个真实的曲率最大点也许很近了，但还不知道有多远，假如这个距离是distance。
2)，可以计算这条曲线的曲率，因为有曲率表达式；但是不可以对表达式求导。在【0,1】区间，曲率变化理论趋势为图2，图3是图1中所有点的曲率，包括A&B，这是因为这个distance还是太大了，漏掉了尖峰，就算有了尖峰还不一定就是曲率k取最大值的点：

图2 L曲线的曲率理论变化趋势，那个尖峰所在区间可能很小

图3 如果分辨率不够就漏掉尖峰了

3），假如二分法去搜索这个点，会让计算量不可接受。


求算法。

mxlzhenzhu

他这个函数不全：http://cda.psych.uiuc.edu/html/4regutools/l_corner.html

还是知识产权问题啊。

看完了，算法有点复杂，自己慢慢搞吧。

mxlzhenzhu

http://www.ncbi.nlm.nih.gov/pubmed/11008433

bazooka529

过来学习了参考了

寂寞的部落

mxlzhenzhu

寂寞的部落 发表于 2013-11-28 16:00

                               
登录/注册后可看大图

gghhjj

1. function [reg_corner,rho,eta,reg_param] = l_curve(U,sm,b,method,L,V)
   
2. %L_CURVE Plot the L-curve and find its "corner".
   
3. %
   
4. % [reg_corner,rho,eta,reg_param] =
   
5. %                  l_curve(U,s,b,method)
   
6. %                  l_curve(U,sm,b,method)  ,  sm = [sigma,mu]
   
7. %                  l_curve(U,s,b,method,L,V)
   
8. %
   
9. % Plots the L-shaped curve of eta, the solution norm || x || or
   
10. % semi-norm || L x ||, as a function of rho, the residual norm
    
11. % || A x - b ||, for the following methods:
    
12. %    method = 'Tikh'  : Tikhonov regularization   (solid line )
    
13. %    method = 'tsvd'  : truncated SVD or GSVD     (o markers  )
    
14. %    method = 'dsvd'  : damped SVD or GSVD        (dotted line)
    
15. %    method = 'mtsvd' : modified TSVD             (x markers  )
    
16. % The corresponding reg. parameters are returned in reg_param.  If no
    
17. % method is specified then 'Tikh' is default.  For other methods use plot_lc.
    
18. %
    
19. % Note that 'Tikh', 'tsvd' and 'dsvd' require either U and s (standard-
    
20. % form regularization) or U and sm (general-form regularization), while
    
21. % 'mtvsd' requires U and s as well as L and V.
    
22. %
    
23. % If any output arguments are specified, then the corner of the L-curve
    
24. % is identified and the corresponding reg. parameter reg_corner is
    
25. % returned.  Use routine l_corner if an upper bound on eta is required.
    
26. 
    
27. % Reference: P. C. Hansen & D. P. O'Leary, "The use of the L-curve in
    
28. % the regularization of discrete ill-posed problems",  SIAM J. Sci.
    
29. % Comput. 14 (1993), pp. 1487-1503.
    
30. 
    
31. % Per Christian Hansen, IMM, July 26, 2007.
    
32. 
    
33. % Set defaults.
    
34. if (nargin==3), method='Tikh'; end  % Tikhonov reg. is default.
    
35. npoints = 200;  % Number of points on the L-curve for Tikh and dsvd.
    
36. smin_ratio = 16*eps;  % Smallest regularization parameter.
    
37. 
    
38. % Initialization.
    
39. [m,n] = size(U); [p,ps] = size(sm);
    
40. if (nargout > 0), locate = 1; else locate = 0; end
    
41. beta = U'*b; beta2 = norm(b)^2 - norm(beta)^2;
    
42. if (ps==1)
    
43.   s = sm; beta = beta(1:p);
    
44. else
    
45.   s = sm(p:-1:1,1)./sm(p:-1:1,2); beta = beta(p:-1:1);
    
46. end
    
47. xi = beta(1:p)./s;
    
48. 
    
49. if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4))
    
50. 
    
51.   eta = zeros(npoints,1); rho = eta; reg_param = eta; s2 = s.^2;
    
52.   reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
    
53.   ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
    
54.   for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
    
55.   for i=1:npoints
    
56.     f = s2./(s2 + reg_param(i)^2);
    
57.     eta(i) = norm(f.*xi);
    
58.     rho(i) = norm((1-f).*beta(1:p));
    
59.   end
    
60.   if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
    
61.   marker = '-'; txt = 'Tikh.';
    
62. 
    
63. elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4))
    
64. 
    
65.   eta = zeros(p,1); rho = eta;
    
66.   eta(1) = abs(xi(1))^2;
    
67.   for k=2:p, eta(k) = eta(k-1) + abs(xi(k))^2; end
    
68.   eta = sqrt(eta);
    
69.   if (m > n)
    
70.     if (beta2 > 0), rho(p) = beta2; else rho(p) = eps^2; end
    
71.   else
    
72.     rho(p) = eps^2;
    
73.   end
    
74.   for k=p-1:-1:1, rho(k) = rho(k+1) + abs(beta(k+1))^2; end
    
75.   rho = sqrt(rho);
    
76.   reg_param = (1:p)'; marker = 'o';
    
77.   if (ps==1)
    
78.     U = U(:,1:p); txt = 'TSVD';
    
79.   else
    
80.     U = U(:,1:p); txt = 'TGSVD';
    
81.   end
    
82. 
    
83. elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4))
    
84. 
    
85.   eta = zeros(npoints,1); rho = eta; reg_param = eta;
    
86.   reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
    
87.   ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
    
88.   for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
    
89.   for i=1:npoints
    
90.     f = s./(s + reg_param(i));
    
91.     eta(i) = norm(f.*xi);
    
92.     rho(i) = norm((1-f).*beta(1:p));
    
93.   end
    
94.   if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
    
95.   marker = ':';
    
96.   if (ps==1), txt = 'DSVD'; else txt = 'DGSVD'; end
    
97. 
    
98. elseif (strncmp(method,'mtsv',4))
    
99. 
    
100.   if (nargin~=6)
     
101.     error('The matrices L and V must also be specified')
     
102.   end
     
103.   [p,n] = size(L); rho = zeros(p,1); eta = rho;
     
104.   [Q,R] = qr(L*V(:,n:-1:n-p),0);
     
105.   for i=1:p
     
106.     k = n-p+i;
     
107.     Lxk = L*V(:,1:k)*xi(1:k);
     
108.     zk = R(1:n-k,1:n-k)\(Q(:,1:n-k)'*Lxk); zk = zk(n-k:-1:1);
     
109.     eta(i) = norm(Q(:,n-k+1:p)'*Lxk);
     
110.     if (i < p)
     
111.       rho(i) = norm(beta(k+1:n) + s(k+1:n).*zk);
     
112.     else
     
113.       rho(i) = eps;
     
114.     end
     
115.   end
     
116.   if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end
     
117.   reg_param = (n-p+1:n)'; txt = 'MTSVD';
     
118.   U = U(:,reg_param); sm = sm(reg_param);
     
119.   marker = 'x'; ps = 2;  % General form regularization.
     
120. 
     
121. else
     
122.   error('Illegal method')
     
123. end
     
124. 
     
125. % Locate the "corner" of the L-curve, if required.
     
126. if (locate)
     
127.   [reg_corner,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,sm,b,method);
     
128. end
     
129. 
     
130. % Make plot.
     
131. plot_lc(rho,eta,marker,ps,reg_param);
     
132. if locate
     
133.   ax = axis;
     
134.   HoldState = ishold; hold on;
     
135.   loglog([min(rho)/100,rho_c],[eta_c,eta_c],':r',...
     
136.          [rho_c,rho_c],[min(eta)/100,eta_c],':r')
     
137.   title(['L-curve, ',txt,' corner at ',num2str(reg_corner)]);
     
138.   axis(ax)
     
139.   if (~HoldState), hold off; end
     
140. end
     

复制代码

1. function [x_k,rho,eta] = mtsvd(U,s,V,b,k,L)
   
2. %MTSVD Modified truncated SVD regularization.
   
3. %
   
4. % [x_k,rho,eta] = mtsvd(U,s,V,b,k,L)
   
5. %
   
6. % Computes the modified TSVD solution:
   
7. %    x_k = V*[ xi_k ] .
   
8. %            [ xi_0 ]
   
9. % Here, xi_k defines the usual TSVD solution
   
10. %    xi_k = inv(diag(s(1:k)))*U(:,1:k)'*b ,
    
11. % and xi_0 is chosen so as to minimize the seminorm || L x_k ||.
    
12. % This leads to choosing xi_0 as follows:
    
13. %    xi_0 = -pinv(L*V(:,k+1:n))*L*V(:,1:k)*xi_k .
    
14. %
    
15. % The truncation parameter must satisfy k > n-p.
    
16. %
    
17. % If k is a vector, then x_k is a matrix such that
    
18. %     x_k = [ x_k(1), x_k(2), ... ] .
    
19. %
    
20. % The solution and residual norms are returned in eta and rho.
    
21. 
    
22. % Reference: P. C. Hansen, T. Sekii & H. Shibahashi, "The modified
    
23. % truncated-SVD method for regularization in general form", SIAM J.
    
24. % Sci. Stat. Comput. 13 (1992), 1142-1150.
    
25. 
    
26. % Per Christian Hansen, IMM, 12/22/95.
    
27. 
    
28. % Initialization.
    
29. m = size(U,1); [p,n] = size(L);
    
30. lk = length(k); kmin = min(k);
    
31. if (kmin<n-p+1 | max(k)>n)
    
32.   error('Illegal truncation parameter k')
    
33. end
    
34. x_k = zeros(n,lk);
    
35. beta = U(:,1:n)'*b; xi = beta./s;
    
36. eta = zeros(lk,1); rho =zeros(lk,1);
    
37. 
    
38. % Compute large enough QR factorization.
    
39. [Q,R] = qr(L*V(:,n:-1:kmin+1),0);
    
40. 
    
41. % Treat each k separately.
    
42. for j=1:lk
    
43.   kj = k(j); xtsvd = V(:,1:kj)*xi(1:kj);
    
44.   if (kj==n)
    
45.     x_k(:,j) = xtsvd;
    
46.   else
    
47.     z = R(1:n-kj,1:n-kj)\(Q(:,1:n-kj)'*(L*xtsvd));
    
48.     z = z(n-kj:-1:1);
    
49.     x_k(:,j) = xtsvd - V(:,kj+1:n)*z;
    
50.   end
    
51.   eta(j) = norm(x_k(:,j));
    
52.   rho(j) = norm(beta(kj+1:n) + s(kj+1:n).*z);
    
53. end
    
54. 
    
55. if (nargout > 1 & m > n)
    
56.   rho = sqrt(rho.^2 + norm(b - U(:,1:n)*beta)^2);
    
57. end

复制代码

1. function [x_k,rho,eta] = tsvd(U,s,V,b,k)
   
2. %TSVD Truncated SVD regularization.
   
3. %
   
4. % [x_k,rho,eta] = tsvd(U,s,V,b,k)
   
5. %
   
6. % Computes the truncated SVD solution
   
7. %    x_k = V(:,1:k)*inv(diag(s(1:k)))*U(:,1:k)'*b .
   
8. % If k is a vector, then x_k is a matrix such that
   
9. %    x_k = [ x_k(1), x_k(2), ... ] .
   
10. %
    
11. % The solution and residual norms are returned in eta and rho.
    
12. 
    
13. % Per Christian Hansen, IMM, 12/21/97.
    
14. 
    
15. % Initialization.
    
16. [n,p] = size(V); lk = length(k);
    
17. if (min(k)<0 | max(k)>p)
    
18.   error('Illegal truncation parameter k')
    
19. end
    
20. x_k = zeros(n,lk);
    
21. eta = zeros(lk,1); rho = zeros(lk,1);
    
22. beta = U(:,1:p)'*b;
    
23. xi = beta./s;
    
24. 
    
25. % Treat each k separately.
    
26. for j=1:lk
    
27.   i = k(j);
    
28.   if (i>0)
    
29.     x_k(:,j) = V(:,1:i)*xi(1:i);
    
30.     eta(j) = norm(xi(1:i));
    
31.     rho(j) = norm(beta(i+1:p));
    
32.   end
    
33. end
    
34. 
    
35. if (nargout > 1 & size(U,1) > p)
    
36.   rho = sqrt(rho.^2 + norm(b - U(:,1:p)*beta)^2);
    
37. end

复制代码