function K = lcrp(p) %LCRP - calculate the Lipschitz constant for the radial projection % % SYNOPSIS: K = lcrp(p) % % Calculates the best Lipschitz constant K for the radial % projection on the space of sequences with the p-norm using % the main result of [1] and described in detail in [2]. % % The value of p should be between 1 and infinity. % % [1] Carlo Franchetti, "The Norm of the Minimum % Projection onto Hyperplanes in L^p[0,1] and the % Radial Constant", Boll. Un. Mat. Ital. 4-B (7), % pp. 803-821, 1990. % [2] J.J. Green, "The Lipschitz constant for the radial % projection on real lp -- implementation notes", % unpublished technical note, 2012. % [3] C. W. Clenshaw, "A note on the summation of Chebyshev % series", Math. Tab. Wash. 9 (1955) 118--120. % [4] Z. Battles & L.N. Trefethen, "An extension of MATLAB to % continuous functions and operators", SIAM J. Sci. Comp., % 25 (2004), 1743--1770. % % Copyright (c) 2011, 2012, 2015, J.J. Green % MIT License if nargin ~= 1 error('exactly one argument required'); end if ~ isscalar(p) error('input must be a scalar'); end if p < 1 error('1 <= p <= infinity required'); end if isinf(p) || (p == 1) K = 2.0; return; end if p == 2 K = 1.0; return; end % if p>2 then coerce p to its conjugate if p > 2 x = 1/(p-1); p = x+1; else x = p-1; end % asymptotic region for x close to zero, see % Section 2 of [2] if x < 1e-9 K = 2 + x * (log(x/8) - 1); return; end % region with x close to 1, see Section 3 of [2] if 1-x < 1e-3 C2 = 0.21961441994532257538; C4 = 0.15444802887153513329; K = polyval([C4 C2 C2 0 1], 1-x); return; end % the Newton-Raphson region 1e-9 < x < 1 - 1e-3, % maximum iterations, typically 2 or 3 are needed. imax = 20; % save a few multiplications from the iteration x2 = x*x; x3 = x2*x; % Initial estimate for maximising u using an order 9 % Chebychev polynomial to approximate the zero-locus % of the derivative, see Section 1 of [2]. These % coefficients generated by Chebfun [4] ucpc = [-2.5595939371497115e-04 +7.0686844943148523e-04 -1.2561575232566727e-03 +1.8047998427177896e-03 -1.5614528129381411e-03 +1.3075424887128736e-03 -1.0205255490737759e-02 +5.6520272597365064e-02 -1.9133303566591867e-01 +2.3504865573520661e-01]; u = chebeval(ucpc, 2*x-1); % initial working values v = u*x; A = v^x; B = v^(1/x); for i = 1:imax xv = x*v; if x > 0.9 % sidestep v^x - v loss of accuracy % as x -> 1 lv = log(v); C = expm1((1/x-1)*lv); D = expm1((x-1)*lv); f = v*( (1+A)*C + x*(1+B)*D ); else f = (1+A)*(B-v) + x*(1+B)*(A-v); end df = (x3 - p*(xv - B)) * A/v ... + (1 - p*(xv - x2*A)) * B/v ... - x*p; % clearly we cannot proceed with zero derivative if abs(df) == 0 error('derivative is zero'); end du = -f/df; u = u + du; % ensure we have not escaped the valid u region if ~ ( (u>0) && (u<1) ) error('iteration has escaped from valid region'); end % update working values (before checking for % convergence, since they are used in the return % value) v = u*x; A = v^x; B = v^(1/x); % stopping criteria % - if |du/u| < eps then u+du = u, nothing more to do % - for second condition see [2] Section 1 if (abs(du) < eps * u) || (abs(f) < x*(1-x)*sqrt(eps)) % typically B is close to zero, so evaluating % (1+B)^x = exp(x * log(1+B)) will lead to loss % of precision, so we use the log1p() function K = ( ((1+A)*exp(x*log1p(B)))^(1/p) ) / (1+v); return; end end error('Newton-Raphson iteration did not converge'); end % Evaluate the Chebychev polynomial with coefficients x % (using the MATLAB convention that the highest-order % coefficient comes first) at the point x, using the % recursive evaluation of Clenshaw [3]. % % The only obscure thing here is the secondary index % variable k, used to index the ring buffer b. Since % MATLAB uses 1-based arrays we need to shift and % unshift the values in the mod call, hence the strange % mod(j-1,3)+1 pattern. % % This function would be easy to vectorise, but this is % not needed for our purposes. function y = chebeval(a, x) N = length(a); if N < 3 error('Chebychev order too small'); end if ~ isscalar(x) error('argument x should be scalar'); end z = 2*x; b = zeros(1, 3); b(1) = a(1); b(2) = z*b(1) + a(2); for j = 1:N-3 k = mod([0:2]+j-1, 3) + 1; b(k(3)) = z*b(k(2)) + a(j+2) - b(k(1)); end k = mod([0:2]+(N-2)-1, 3) + 1; y = x*b(k(2)) + a(N) - b(k(1)); end