As a follow-up on the previous post titled Projection on an epigraph, we here discuss how we can project on the epigraph of the squared norm function.

Define with and take .

Let . Following the previous post, the optimality conditions are

Then, the first optimality condition becomes

Then,

where . This third-order polynomial equation can be solved numerically, it has a root which satisfies the second optimality condition above. Once has been determined, can be computed.

The above can be trivially extended to the case where , with $\alpha_i>0$. However, the projection onto the epigraph of where is a symmetric positive (semi)definite matrix is more complex and we cannot derive explicit formulas in the general case.

Here is a simple MATLAB function to compute the projection onto the epigraph of the squared Euclidean norm efficiently:

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

### Like this:

Like Loading...