[guided]Fix one function $f\in\mathcal F$. The purpose of this step is to get a deviation estimate for that single function before taking a union over the class. For each $i\in\{1,\dots,n\}$, define
\begin{align*}
X_{i,f}:\Omega\to\mathbb R,\quad X_{i,f}(\omega)=f(Z_i(\omega))-\mathbb E[f(Z_i)].
\end{align*}
These are the centered observations associated with $f$. Define their empirical average by
\begin{align*}
M_f:\Omega\to\mathbb R,\quad M_f(\omega)=\frac{1}{n}\sum_{i=1}^n X_{i,f}(\omega).
\end{align*}
We want to bound $\mathbb P(|M_f|>t)$ for a number $t>0$.
First bound the upper tail. Let $\lambda>0$. Since the exponential function is increasing,
\begin{align*}
\{M_f\ge t\}
=
\left\{\exp(\lambda nM_f)\ge \exp(\lambda nt)\right\}.
\end{align*}
For any nonnegative random variable $Y$ and number $a>0$, the inequality $a\,\mathbb 1_{\{Y\ge a\}}\le Y$ implies $\mathbb P(Y\ge a)\le \mathbb E[Y]/a$. Applying this with $Y=\exp(\lambda nM_f)$ and $a=\exp(\lambda nt)$ gives
\begin{align*}
\mathbb P(M_f\ge t)
&\le
\exp(-\lambda nt)\mathbb E[\exp(\lambda nM_f)].
\end{align*}
Now compute the exponential moment. Since
\begin{align*}
nM_f=\sum_{i=1}^n X_{i,f},
\end{align*}
and the random variables $X_{1,f},\dots,X_{n,f}$ are independent for this fixed $f$, the expectation factors:
\begin{align*}
\mathbb E[\exp(\lambda nM_f)]=\mathbb E\left[\exp\left(\lambda\sum_{i=1}^n X_{i,f}\right)\right].
\end{align*}
Since the exponential of a sum is the product of exponentials,
\begin{align*}
\mathbb E\left[\exp\left(\lambda\sum_{i=1}^n X_{i,f}\right)\right]=\mathbb E\left[\prod_{i=1}^n \exp(\lambda X_{i,f})\right].
\end{align*}
Independence of $X_{1,f},\dots,X_{n,f}$ gives
\begin{align*}
\mathbb E\left[\prod_{i=1}^n \exp(\lambda X_{i,f})\right]=\prod_{i=1}^n \mathbb E[\exp(\lambda X_{i,f})].
\end{align*}
The sub-Gaussian hypothesis gives, for every $i$,
\begin{align*}
\mathbb E[\exp(\lambda X_{i,f})]\le \exp\left(\frac{\lambda^2\sigma^2}{2}\right).
\end{align*}
Therefore
\begin{align*}
\mathbb E[\exp(\lambda nM_f)]\le \prod_{i=1}^n \exp\left(\frac{\lambda^2\sigma^2}{2}\right)=\exp\left(\frac{n\lambda^2\sigma^2}{2}\right).
\end{align*}
Substituting this into the upper-tail estimate yields
\begin{align*}
\mathbb P(M_f\ge t)
\le
\exp\left(-\lambda nt+\frac{n\lambda^2\sigma^2}{2}\right).
\end{align*}
The exponent is minimized over $\lambda>0$ by choosing $\lambda=t/\sigma^2$, which is positive because $t>0$ and $\sigma>0$. With this choice,
\begin{align*}
-\lambda nt+\frac{n\lambda^2\sigma^2}{2}=-\frac{nt^2}{\sigma^2}+\frac{nt^2}{2\sigma^2}=-\frac{nt^2}{2\sigma^2}.
\end{align*}
Hence
\begin{align*}
\mathbb P(M_f\ge t)\le \exp\left(-\frac{nt^2}{2\sigma^2}\right).
\end{align*}
The lower tail is the same exponential argument applied to $-M_f$. Equivalently, use the sub-Gaussian moment bound with $-\lambda$ in place of $\lambda$:
\begin{align*}
\mathbb E[\exp(-\lambda X_{i,f})]\le \exp\left(\frac{\lambda^2\sigma^2}{2}\right).
\end{align*}
Repeating the preceding computation gives
\begin{align*}
\mathbb P(M_f\le -t)\le \exp\left(-\frac{nt^2}{2\sigma^2}\right).
\end{align*}
Combining the upper and lower tails by subadditivity of probability,
\begin{align*}
\mathbb P(|M_f|>t)=\mathbb P\bigl(\{M_f>t\}\cup\{M_f<-t\}\bigr).
\end{align*}
By subadditivity of probability,
\begin{align*}
\mathbb P\bigl(\{M_f>t\}\cup\{M_f<-t\}\bigr)\le \mathbb P(M_f>t)+\mathbb P(M_f<-t)\le 2\exp\left(-\frac{nt^2}{2\sigma^2}\right).
\end{align*}
This is the fixed-function concentration estimate.[/guided]