<s>[INST] How to properly find supremum of a function $f(x,y,z)$ on a cube $[0,1]^3$? Solving an applied problem I was faced with the need to find supremum of the following function 
$$f(x,y,z)=\frac{(x-xyz)(y-xyz)(z-xyz)}{(1-xyz)^3}$$
where $f\colon\ [0,1]^3\backslash\{(1,1,1)\} \to\mathbb{R}$.
I used Wolfram Mathematica and it resulted with "there is no maximum in the region in which the objective function is defined" and that the supremum is $\frac{8}{27}$ when $x,y,z\to 1$. I was able to verify some of it:
$$\lim\limits_{x,y,z\to 1}f(x,y,z) 
= 
\lim\limits_{x\to 1}f(x,x,x)
=
\lim\limits_{x\to 1}\frac{(x-x^3)^3}{(1-x^3)^3}
=
\lim\limits_{x\to 1}\frac{x^3(x+1)^3}{(x^2+x+1)^3}
=
\frac{8}{27}
$$ 
though I'm not sure how to justify the first equality. Can I neglect direction of derivative in this case?
And secondly, I don't know how to show that $f(x,y,z)< \frac{8}{27}$ for all $x,y,z\in[0,1]^3\backslash\{(1,1,1)\}$. When I tried to do that straightforwardly I drowned in computations. 
 [/INST] tl; dr: The supremum of $f$ can be found by restricting to the diagonal of the cube and approaching $(1, 1, 1)$.

Partition the cube $C$ into level sets of $xyz$. If $0 < c < 1$, the level set
$$
L_{c} = \{(x, y, z) \in C : f(x, y, z) = c\}
$$
is compact, so $f$ achieves absolute extrema in $L_{c}$. (In fact, $L_{c}$ is homeomorphic to a closed disk, with three hyperbola arcs as boundary.)
To establish the bound
$$
f(x, y, z) \leq f(\sqrt[3]{xyz}, \sqrt[3]{xyz}, \sqrt[3]{xyz})
  = f(\sqrt[3]{c}, \sqrt[3]{c}, \sqrt[3]{c})\quad\text{on $L_{c}$,}
\tag{1}
$$
it suffices to show that the maximum of $f$ subject to $xyz = c$ occurs where $x = y = z$.
Assume from now on that $0 < c < 1$. Introduce the function $g(u) = u - \frac{1}{u}$ with domain $0 < u < 1$, and note that $ug'(u) = u + \frac{1}{u}$ is monotone for $0 < u < 1$.
On the surface $xyz = c$, we have $xy = \frac{c}{z}$, $xz = \frac{c}{y}$, and $yz = \frac{c}{x}$, so
\begin{align*}
f(x, y, z)
  &= \frac{(x - c)(y - c)(z - c)}{(1- c)^{3}} \\
  &= \frac{c - c(yz + xz + xy) + c^{2}(x + y + z) - c^{3}}{(1- c)^{3}} \\
  &= \frac{c - c^{2}(\frac{1}{x} + \frac{1}{y} + \frac{1}{z}) + c^{2}(x + y + z) - c^{3}}{(1- c)^{3}} \\
  &= \frac{c + c^{2}\bigl(g(x) + g(y) + g(z)\bigr) - c^{3}}{(1- c)^{3}}. \\
\end{align*}
Since the denominator is positive, this is maximized when $g(x) + g(y) + g(z)$ is maximized subject to $xyz = c$.
Lagrange multiplers gives $g'(x) = \lambda yz$, or $xg'(x) = c\lambda$, and similarly for $y$ and $z$. That is, the only critical points of $f(x, y, z)$ on the surface $xyz = c$ occur where
$$
xg'(x) = yg'(y) = zg'(z).
$$
Since $0 < x, y, z < 1$ on the cube and $u \mapsto ug'(u)$ is monotone (hence injective) in $(0, 1)$, the preceding implies $x = y = z$. That is, the only critical point of $f$ in the interior of $L_{c}$ (and hence the only possible interior extremum) is $(\sqrt[3]{xyz}, \sqrt[3]{xyz}, \sqrt[3]{xyz}) = (\sqrt[3]{c}, \sqrt[3]{c}, \sqrt[3]{c})$.
To conclude (1) holds, it remains to show that the boundary values of $f$ do not exceed the right-hand side. Because $f$ is symmetric, we may as well assume $x = 1$ and $yz = c$, so that
$$
f(x, y, z)
  = \frac{(y - c)(z - c)}{(1- c)^{2}}
  = \frac{c - c(y + z) + c^{2}}{(1 - c)^{2}}.
$$
This is maximized when $y + z$ is minimized, which is easily seen to occur where $y = z = \sqrt{c}$. Since
$$
f(1, \sqrt{c}, \sqrt{c}) < f(\sqrt[3]{c}, \sqrt[3]{c}, \sqrt[3]{c}),
$$
(a straightforward, if slightly tedious, one-variable calculation), (1) holds.
</s>