Projection onto the epigraph of the squared norm

epipr_sqnorm.m

function [x_, z_, details] = epipr_sqnorm(x,z)
if (x'*x <= z)
    x_ = x; z_ = z;
    return
end
theta = 1 - 2 * z;
[r, status] = cubic_roots(theta, x);
details.status = status;
% Pick the right root
for i=1:length(r),   
    x_ = x/(1 + 2*(r(i) - z));
    if abs(norm(x_)^2-r(i)) < 1e-6,
        z_ = r(i);
        break;
    end
end
% Refine the root
[z_, iter, err] = newton_solver(x,theta,z_,5,1e-10);
x_ = x/(1 + 2*(z_ - z));
details.newton_iter = iter;
details.newton_err  = err;

function [r, status] = cubic_roots(theta, x)
b=4*theta; c=theta^2; d=-x'*x;
D  = 72*b*c*d - 4*b^3*d +b^2*c^2 - 16*c^3 - 432*d^2;
D0 = b^2 - 12*c;
status.D = D; status.D0 = D0;
if abs(D)<1e-14,
    if abs(D0)<1e-14,
        % one triple root (we cannot be here!)
        r = -b/12;
        status.msg = 'one triple root';
    else
        % a double root and a single one
        r = zeros(2,1);
        r(1) = (16*b*c - 144*d - b^3)/(4*D0); % single
        r(2) = (36*d - b*c)/(2*D0); % double (cannot be)
        status.msg = 'double plus single';
    end
    return;
end
r = roots([4 b c d]); % eigenvalues of matrix

function [zsol, i, err] = ...
   newton_solver(x,theta,z0,maxit,tolx)
zsol = z0; zsol_prev = z0-1;
for i=1:maxit,
  zsol = zsol - ...
    (4*zsol ^3 +  4*theta*zsol^2 + theta^2 * zsol -x'*x)...
    /(12*zsol^2 + 8*theta*zsol + theta^2);
  err = abs(zsol-zsol_prev);
  if err<tolx,
    return;
  end
  zsol_prev = zsol;
end