The line Aubin says to "develop" is the inner product of the $(0, 3)$ tensor $T_{\nu\alpha\beta} = \nabla_\nu \nabla_\alpha \psi \nabla_\beta \psi - \nabla_\nu \nabla_\beta \psi \nabla_\alpha \psi$ with itself, i.e. $|T|^2_{g}$. That's why it's positive.
By "develop" Aubin means "distribute out the multiplication", which lands us at
\begin{align}
|T|^2_{g}
&=
g^{\nu \mu}g^{\alpha \lambda}g^{\beta \gamma}
\nabla_\nu \nabla_\alpha \psi \nabla_\beta \psi \nabla_{\mu}\nabla_\lambda\psi\nabla_{\gamma}\psi \\
&-
g^{\nu \mu}g^{\alpha \lambda}g^{\beta \gamma}
\nabla_\nu \nabla_\alpha \psi \nabla_\beta \psi
\nabla_\mu\nabla_\gamma \psi \nabla_\lambda \psi \\
&-
g^{\nu \mu}g^{\alpha \lambda}g^{\beta \gamma}
\nabla_\nu \nabla_{\beta} \psi \nabla_\alpha {\psi}
\nabla_{\mu}\nabla_\lambda\psi\nabla_{\gamma}\psi \\
&+
g^{\nu \mu}g^{\alpha \lambda}g^{\beta \gamma}
\nabla_\nu \nabla_{\beta} \psi \nabla_\alpha {\psi}
\nabla_\mu\nabla_\gamma \psi \nabla_\lambda \psi.
\end{align}
We can recognize the first and last lines as $|\nabla \nabla \psi|^2|\nabla \psi|^2$.
For the other lines, it's easier to do the calculation backwards: note
\begin{align}
|\nabla |\nabla \psi|^2|^2
&= g^{ab}
\nabla_a(g^{cd}\nabla_c\psi\nabla_d\psi)
\nabla_b(g^{ef}\nabla_e\psi\nabla_f\psi)\\
\end{align}
since I can switch $c,d$ in the inverse-metric we get the same term when the $\nabla_a$ hits either $\nabla_c$ or $\nabla_d$. Similarly with $\nabla_b$. So
\begin{align}
|\nabla |\nabla \psi|^2|^2
&= 4g^{ab}g^{cd}g^{ef}
(\nabla_a\nabla_c\psi\nabla_d\psi)
(\nabla_b\nabla_e\psi\nabla_f\psi)\\
\end{align}
Now, you can match the terms above with the two terms which come with a minus sign in our big calculation of $|T|^2_g$. We find,
$$|T|^2_g = 2|\nabla \nabla \psi|^2 |\nabla \psi|^2 - (1/2)|\nabla |\nabla \psi|^2|^2.$$
This proves the claim at the end.
To use the claim, let $f = |\nabla \psi|$ and note
$$|T|^2_g = 2|\nabla \nabla \psi|^2 f^2 - (1/2)|\nabla f^2|^2$$
Now do the derivative, divided by $f^2$, and find what you want.
If you're reading this inequality for the first time, it may be useful to note the exact constant of $1$ in the inequality is sometimes important in geometric analysis for elliptic pde. It comes up when calculating the laplacian of norms of $u$ or norms of derivatives of $u$, assuming you know something about the laplacian of $u$. Say $u$ is a tensor, if we calculate
$\Delta |u|^2$ we get a positive term $2|\nabla u|^2$. If you then use that to calculate $\Delta |u|$ you pick up a negative term (since $|u| = \sqrt{|u|^2}$ and $x \mapsto \sqrt{x}$ is concave) which is always beat by the positive term $2|\nabla u|^2$, using the inequality.