Edit: I realized that the explanation of the former second step of proof below was a little bit obscure since, while entirely correct, did not clarify enough why the choice of integrability exponent is not done by guessing. Therefore I decided to substitute it by a similar but more direct procedure and put the former step 2 in the notes for a brief proof of their equivalence.
I have not seen the proof of Young's inequality you allude to: however, the answer to your question i.e. whether it is possible to prove it in another way, remaining at the same level of knowledge, is yes. The proof given below is inspired by and follows the one in the nice monograph [1], pp. 26-27 and it is based on the standard and generalized (i.e. involving three or more functions) Hölder's inequalities and by a judicious choice of the integrability exponents associated to two or three factors expressing $| f(y)g(x-y) |$: the introduction of an auxiliary function $h$ is not required.
Young's Inequality. Let $p,q,r\in\Bbb R$ be such that
$$
1\le p\le q\le +\infty, \quad 1+\frac{1}{r}=\frac{1}{p}+\frac{1}{q},\label{1}\tag{1}
$$
and let $f\in L^p(\mathbf{R}^d)$ and $g\in L^q(\mathbf{R}^d)$: if
$$
f\ast g(x)=\int\limits_{\mathbf{R}^d} f(y)g(x-y)\mathrm{d}y,
$$
then
$$
\Vert f\ast g\Vert_{r}\le \Vert f\Vert_{p}\Vert g\Vert_q\label{2}\tag{2}
$$
Comment. The strategy of the proof goes as follows: first, in every range of values of $p, q, r$ defined by conditions \eqref{1}, we will express $|f(y)g(x-y)|$ as the product of three factors
$$
| f(y)g(x-y) |=\big(|f(y)|^p|g(x-y)|^q\big)^\frac{1}{s_1}|g(x-y)|^{\frac{q}{s_2}} |f(y)|^{\frac{p}{s_3}}.\label{step1}\tag{Step 1}
$$
Since we want to estimate the the $L^r$ norm of the convolution, we assume $s_1=r$.
Now \ref{step1} implies that the coefficients $s_1, s_2, s_3$ must satisfy the following conditions
$$
p\left(\frac{1}{s_1}+\frac{1}{s_3}\right)=1\quad q\left(\frac{1}{s_1}+\frac{1}{s_2}\right)=1,\label{c1}\tag{C1}
$$
We thus have a non-homogeneous linear system in the $s_i^{-1}$ variables, $i=1,2,3$ which is uniquely solvable, provided $pq\neq0$, and the second step consist in solving it for the unknown exponents: explicitly
$$
\begin{pmatrix}
1 & 0 & 0\\
q & q & 0\\
p & 0 & p
\end{pmatrix}
\begin{pmatrix}
s_1^{-1}\\
s_2^{-1}\\
s_3^{-1}\\
\end{pmatrix}=
\begin{pmatrix}
{\frac{1}{r}}\\
{1}\\
{1}\\
\end{pmatrix}\iff
\begin{pmatrix}
s_1^{-1}\\
s_2^{-1}\\
s_3^{-1}\\
\end{pmatrix}=
\begin{pmatrix}
{\frac{1}{r}}\\
{\frac{1-\frac{q}{r}}{q}}\\
{\frac{1-\frac{p}{r}}{p}}\\
\end{pmatrix}\label{step2}\tag{Step 2}
$$
The third and final step is to estimate the $L^r$ norm of the convolution $f\ast g$ by applying to equation \eqref{step1} one of the various forms of Hölder's inequality. This of course can be done since it is easily verified that
$$
\frac{1}{s_1} + \frac{1}{s_2} + \frac{1}{s_3}=1.\label{c2}\tag{C2}
$$
Proof. If $r=\infty$, then \eqref{2} is a direct consequence of the standard Holder's inequality, since
$$
\frac{1}{p}+\frac{1}{q}=1.
$$
Assuming $r<+\infty$, \eqref{2} must be verified for the three ranges defined by conditions \eqref{1}, i.e.
- $1<p<r$ and $1<q< r$;
- $p=1<q=r$;
- $p=r$ and $q=1$.
- Case 1: this is the most general case. From \ref{step2} we have
$$
\begin{cases}
s_1=r\\
\\
s_2=\dfrac{q}{1-\frac{q}{r}}\\
s_3=\dfrac{p}{1-\frac{p}{r}}
\end{cases},
$$
and thus equation \eqref{step1} becomes
$$
| f(y)g(x-y) |=\big(|f(y)|^p|g(x-y)|^q\big)^\frac{1}{r}|g(x-y)|^{1-\frac{q}{r}} |f(y)|^{1-\frac{p}{r}}.\label{3}\tag{3}
$$
Estimating the convolution $f\ast g$ by using \eqref{3} and the generalized Hölder inequality gives
$$
|f\ast g(x)|\le \bigg(\int\limits_{\mathbf{R}^d} |f(y)|^p|g(x-y)|^q\mathrm{d}y\bigg)^{\!\frac{1}{r}}\Vert f\Vert_{p}^{1-\frac{p}{r}} \Vert g\Vert_q^{1-\frac{q}{r}}\label{4}\tag{4}
$$
and applying the generalized Hölder inequality to \eqref{4} finally gives
$$
\begin{split}
\Vert f\ast g\Vert_r &\le \Vert f\Vert_{p}^{1-\frac{p}{r}} \Vert g\Vert_r^{1-\frac{q}{r}}\bigg(\int\limits_{\mathbf{R}^d}\mathrm{d}x \int\limits_{\mathbf{R}^d} |f(y)|^p|g(x-y)|^q\mathrm{d}y\bigg)^{\frac{1}{r}}\\
& = \Vert f\Vert_{p}^{1-\frac{p}{r}} \Vert g\Vert_q^{1-\frac{q}{r}} \bigg(\int\limits_{\mathbf{R}^d}|f(y)|^p \mathrm{d}y \int\limits_{\mathbf{R}^d} |g(x)|^q\mathrm{d}x\bigg)^{\!\frac{1}{r}}\\
& = \Vert f\Vert_{p}\Vert g\Vert_r
\end{split}
$$
- Case 2 and Case 3: in these cases, the right side of equation \eqref{3} reduces to the product of two terms and inequality \eqref{2} is obtained by means of the standard Hölder inequality. Explicitly,
$$
| f(y)g(x-y)| =
\begin{cases}
\big(|f(y)||g(x-y)|^q\big)^\frac{1}{r}|f(y)|^{1-\frac{1}{r}}&\text{ in case 2}\\
\big(|f(y)|^p|g(x-y)|\big)^\frac{1}{r}|g(x-y)|^{1-\frac{1}{r}}&\text{ in case 3}
\end{cases}.\qquad\blacksquare
$$
Final notes
- This proof is entirely elementary and does not require the "guessing" of the coefficients $s_1, s_2, s_3$, which are instead well defined and calculable.
- Why the second step above is entirely equivalent to the formerly proposed one? Because the non-homogeneous linear system obtained by considering directly, without assuming a priori $s_1=r$, the conditions implied by equation \eqref{3} and by the necessity of using Hölder's inequality i.e. \eqref{c1} and \eqref{c2}, is perfectly equivalent to \eqref{step2}. To see this, is sufficient to write it down and solve it
$$
\begin{pmatrix}
p & 0 & p\\
q & q & 0\\
1 & 1 & 1
\end{pmatrix}
\begin{pmatrix}
s_1^{-1}\\
s_2^{-1}\\
s_3^{-1}\\
\end{pmatrix}=
\begin{pmatrix}
{1}\\
{1}\\
{1}\\
\end{pmatrix}\iff
\begin{pmatrix}
s_1^{-1}\\
s_2^{-1}\\
s_3^{-1}\\
\end{pmatrix}=
\begin{pmatrix}
{\frac{1}{p}+\frac{1}{p}-1}\\
{1-\frac{1}{p}}\\
{1-\frac{1}{q}}\\
\end{pmatrix},
$$
and then use \eqref{1} to express $s_1, s_2$ and $s_3$ respectively as functions of $r$, $r$ and $q$, $r$ and $p$. Finally, it is wort to note that, in the approach above, the truth of \eqref{c2} is a consequence of the implicit use of \eqref{1}.
- As a final remark, let me say that Besov, Il'in and Nikol’skiĭ prove \eqref{2} ([1], p. 27-28) first for $d=1$ (without showing explicitly \eqref{step1} and \eqref{step2}) and then for vector exponents and $d\ge 2$, i.e. $\mathbf{p}=(p_1,\ldots,p_d)$, $\mathbf{q}=(q_1,\ldots,q_d)$ and $\mathbf{r}=(r_1,\ldots,r_d)$ where each of their $i$-th component satisfies relation \eqref{1}: the result is used in the development of the theory of anysotropic function (Sobolev and Besov) spaces.
Bibliography
[1] Oleg V. Besov, Valentin P. Il’in, Sergei M. Nikol’skiĭ (1978), Integral representations of functions and imbedding theorems. Vol. I, Ed. by Mitchell H. Taibleson. Translation from the Russian. (English) Scripta Series in Mathematics. Washington, D.C.: V. H. Winston & Sons. New York-Toronto-London: John Wiley & Sons, ISBN: 0-470-26540-X, pp. VIII+345, MR0519341, Zbl 0392.46022.
Best Answer
This inequality cannot hold in general. Indeed, \begin{equation*} |\nabla(f^{\frac{p+q-2}{p}})|^p=k^pf^{q-2}|\nabla f|^p, \end{equation*} where \begin{equation*} k:=\frac{p+q-2}p=1+\frac{q-2}p\ge1. \end{equation*} So, at all points where $g>0$, $|\nabla g|>0$, and $|\nabla f|>0$, we can rewrite the inequality in question as \begin{equation*} r-1\le C(s,q)(r^p-k^p) \tag{1}\label{1} \end{equation*} where $r:=a/b$, $a:=f/g$, and $b:=|\nabla f|/|\nabla g|$.
In general, $r$ can take any nonnegative real value.
Letting now $r\to\infty$ in \eqref{1}, we get $C(s,q)>0$. Letting then $r=1$, we get a contradiction: $0\le C(s,q)(1-k)<0$ if $k>1$, that is, if $q>2$.
If, finally, $q=2$, then $k=1$ and the only value of $C(s,q)$ such that \eqref{1} holds for all real $r\ge0$ is $1/p$ -- so that $C(s,2)$ must depend, not on $s$, but on $p$.