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}


\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;
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);
% 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';
    % 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';
r = roots([4 b c d]); % eigenvalues of matrix

function [zsol, i, err] = ...
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,
  zsol_prev = zsol;



Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s


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 is the mathematician of the village of Asterix

%d bloggers like this: