# Projection on the epigraph of the squared Euclidean norm

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 $f:\mathbb{R}^n\to\mathbb{R}:x\mapsto \|x\|^2$ with $\nabla f(x) = 2x$ and take $(z,\tau)\notin \mathrm{epi} f$.

Let $\Pi_f(z,\tau) = (\bar{x},\bar{s})$. Following the previous post, the optimality conditions are

\begin{aligned} z - \bar{x} &= 2(\|\bar{x}\|^2 - \tau)\bar{x},\\ \bar{s} &= \|\bar{x}\|^2. \end{aligned}

Then, the first optimality condition becomes

\begin{aligned} & z - \bar{x} = 2(\bar{s} - \tau)\bar{x}\notag \\ \Leftrightarrow\ & \bar{x} = \frac{z}{1+2\bar{s}-2\tau}\label{eq:sqnorm:barx} \end{aligned}

Then,

\begin{aligned} &\bar{s} = \|\bar{x}\|^2 = \left\| \frac{z}{1+\bar{s}-\tau} \right\|^2 = \frac{\|z\|^2}{(1+2\bar{s}-2\tau)^2} \notag \\ \Leftrightarrow\ & \bar{s}(1+2\bar{s}-2\tau)^2 = \|z\|^2\\ \Leftrightarrow\ & 4\bar{s}^3 + 4\theta \bar{s}^2 + \theta^2 \bar{s} - \|z\|^2 = 0, \end{aligned}

where $\theta = 1-2\tau$. This third-order polynomial equation can be solved numerically, it has a root $\bar{s}=\bar{s}(\theta, \|z\|^2)$ which satisfies the second optimality condition above. Once $\bar{s}$ has been determined, $\bar{x}$ can be computed.

The above can be trivially extended to the case where $f(x) = \sum_{i=1}^{n}\alpha_i x_i^2$, with $\alpha_i>0$. However, the projection onto the epigraph of $f(x) = x'Px$ where $P$ 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


mathbabe

Exploring and venting about quantitative issues

Look at the corners!

The math blog of Dmitry Ostrovsky

The Unapologetic Mathematician

Mathematics for the interested outsider

Almost Sure

A random mathematical blog

Mathematix

Mathematix is the mathematician of the village of Asterix