Given a nonempty closed set $C\subseteq\mathbb{R}^n$, let $S_C\subseteq\mathbb{R}^n$ denote the set of points for which the closest point in $C$ is not unique. Suppose $C$ satisfies each of the following:
- $S_C$ is nonempty (i.e., $C$ is not convex),
- there exists $\rho>0$ such that $C^\rho\cap S_C=\emptyset$, and
- the closest point mapping $c\colon(\mathbb{R}^n\setminus S_C)\to C$ is continuous.
(Every nonempty closed nonconvex subset $C\subseteq\mathbb{R}^n$ with smooth boundary that I can think of satisfies these property.)
Let $r_C<\infty$ denote the supremum of such $\rho$. Below, we show that for every $\epsilon\in(0,r_C)$, it holds that $u_C$ is not Lipschitz on $\mathbb{R}^n\setminus (C^\epsilon\cup Z)$ for any null set $Z$.
Select $p\in\mathbb{R}^n$ of distance $r_C$ from $C$ for which there exist $x,y\in C$ with $x\neq y$ such that
$$
\|x-p\|=\|y-p\|=r.
$$
Then for each $t\in(0,1)$, the points $x(t):=tx+(1-t)p$ and $y(t):=ty+(1-t)p$ are of distance less than $r_C$ from $C$, and so $c(x(t))=x$ and $c(y(t))=y$. Furthermore,
$$
u_C(x(t))=\frac{x(t)-c(x(t))}{\|x(t)-c(x(t))\|}=\frac{x(t)-x}{\|x(t)-x\|}=\frac{p-x}{\|p-x\|},
$$
and similarly
$$
u_C(y(t))=\frac{p-y}{\|p-y\|}.
$$
Notably, for small $t$, $x(t)$ and $y(t)$ are arbitrarily close to each other, while $u_C(x(t))$ and $u_C(y(t))$ remain a fixed distance apart.
Given $\epsilon\in(0,r_C)$, a null set $Z$, and a constant $L>0$, we now show that $u_C$ is not $L$-Lipschitz on $\mathbb{R}^n\setminus (C^\epsilon\cup Z)$. Fix $\delta>0$ to be selected later. Select $t_0\in(0,1/2)$ small enough so that $x(t_0)$ and $y(t_0)$ avoid $C^\epsilon$ and $\|x(t_0)-y(t_0)\|<\delta$. Select a continuous path $q\colon[0,1]\to \mathbb{R}^n\setminus S_C$ such that $q(0)=x(t_0)$, $q(1)=y(t_0)$, and $q(s)\not\in Z$ for almost every $s\in(0,1)$. By the continuity of $c$, there exist $s_0,s_1\in[0,1]$ such that $q(s_0),q(s_1)\not\in Z$ and
$$
\max\Big\{\|q(s_0)-x(t_0)\|,\|q(s_1)-y(t_0)\|,\|c(q(s_0))-x\|,\|c(q(s_1))-y\|\Big\}<\delta.
$$
As we will see, $q(s_0)$ and $q(s_1)$ witness that $u_C$ is not $L$-Lipschitz on $\mathbb{R}^n\setminus (C^\epsilon\cup Z)$.
After some straightforward manipulations, we have
\begin{align*}
\|u_C(q(s_0))-u_C(x(t_0))\|
&=\bigg\|\frac{q(s_0)-c(q(s_0))}{\|q(s_0)-c(q(s_0))\|}-\frac{x(t_0)-x}{\|x(t_0)-x\|}\bigg\|\\
&\leq\frac{4\delta}{\|x(t_0)-x\|}=\frac{4\delta}{(1-t_0)r_C}\leq\frac{8\delta}{r_C},
\end{align*}
and similarly
$$
\|u_C(q(s_1))-u_C(y(t_0))\|\leq\frac{8\delta}{r_C}.
$$
Then
$$
\|u_C(q(s_0))-u_C(q(s_1))\|
\geq\|u_C(x(t_0))-u_C(y(t_0))\|-\frac{16\delta}{r_C}
=\frac{\|x-y\|-16\delta}{r_C},
$$
but
$$
\|q(s_0)-q(s_1)\|
\leq\|x(t_0)-y(t_0)\|+2\delta
\leq 3\delta.
$$
Selecting $\delta<\min\{\frac{1}{32},\frac{1}{6Lr_c}\}\cdot\|x-y\|$ then gives
$$
\|u_C(q(s_0))-u_C(q(s_1))\|>L\|q(s_0)-q(s_1)\|,
$$
as claimed.
Seems like in fact convex sets are the only ones for which $p_C$ is continuous. To prove this, we can begin by noticing that for a set $C$ with the property, $p_C$ can only take one point sets as values. If not, let $x,a_1,a_2$ be different points with $a_1,a_2\in p_C(x)$. Then any point $y$ in the open segment from $x$ to $a_i$ has $p_C(y)=\{a_i\}$, so $p_C$ cannot be continuous at $x$.
So, if for a compact $C$, $p_C$ is continuous, then it is single valued. So this shows that $C$ has to be convex.
Edit: From Dustin G. Mixon's comment above, it seems like the nearest point mapping is well studied. A reference that if a compact set $C$ is not convex, then the nearest point mapping is not single valued would be welcome.
Best Answer
EDIT: The OP changed the question while I was writing this answer. In this answer, we take the domain of $u$ to be $\{x\in\mathbb{R}^n:d(x)\geq\epsilon\}$ for some fixed $\epsilon>0$. This is a different use of $\epsilon$ than in OP's current problem formulation.
Let $B_\epsilon$ denote the open ball centered at the origin of radius $\epsilon$. Below, we show how to factor $u$ as the composition of a $2$-Lipschitz map $g\colon\mathbb{R}^n\to\mathbb{R}^n$ with a $(1/\epsilon)$-Lipschitz map $f\colon(\mathbb{R}^n\setminus B_\epsilon)\to S^{n-1}$.
To analyze each map, we will use the fact that projection onto a nonempty closed convex set is contractive.
First, define $g$ by $g(x):=x-p(x)$, where $p$ denotes projection onto $P$. Then $$ \|g(x)-g(y)\|=\|(x-p(x))-(y-p(y))\|\leq\|x-y\|+\|p(x)-p(y)\|\leq2\|x-y\|. $$ Next, define $f$ by $f(x):=x/\|x\|$. Since the domain of $f$ is the complement of $B_\epsilon$, we can instead write $f(x)=q(x)/\epsilon$, where $q$ denotes projection onto the closure of $B_\epsilon$. Then $$ \|f(x)-f(y)\|=\|q(x)/\epsilon-q(y)/\epsilon\|=(1/\epsilon)\|q(x)-q(y)\|\leq(1/\epsilon)\|x-y\|. $$ Combining these estimates, we have that $u=f\circ g$ is $(2/\epsilon)$-Lipschitz.