[proofplan]
We prove an upper-tail estimate for the largest eigenvalue of $S:=\sum_{i=1}^nY_i$ by applying the matrix [Laplace transform](/page/Laplace%20Transform) method to $\operatorname{tr}\exp(\theta S)$. The scalar Bernstein exponential bound gives a Loewner-order estimate for each matrix moment generating function, while the noncommutative trace moment generating function estimate combines the independent summands inside a single trace. Optimising the resulting one-parameter exponential bound gives the Bernstein denominator $\sigma^2+Lt/3$. Applying the same argument to $-S$ gives the lower tail, and the operator norm bound follows by the union bound.
[/proofplan]
[step:Reduce to the non-degenerate case and introduce the variance matrix]
Let $(\Omega,\mathcal F,\mathbb P)$ denote the probability space on which the matrices $Y_1,\dots,Y_n$ are defined, and let $\mathbb E$ denote expectation with respect to $\mathbb P$. Let $\mathbb S^d\subset \mathbb R^{d\times d}$ denote the [vector space](/page/Vector%20Space) of real symmetric $d\times d$ matrices, and let $I\in\mathbb S^d$ denote the $d\times d$ identity matrix. For a matrix $A\in\mathbb S^d$, let $\lambda_{\max}(A)$ denote its largest eigenvalue, let $\operatorname{tr}(A)$ denote its trace, and let $\exp(A)\in\mathbb S^d$ denote the matrix exponential defined by the absolutely convergent [power series](/page/Power%20Series)
\begin{align*}
\exp(A)=\sum_{m=0}^{\infty}\frac{A^m}{m!}.
\end{align*}
Define the random matrix $S:\Omega\to \mathbb S^d$ by
\begin{align*}
S(\omega)=\sum_{i=1}^nY_i(\omega).
\end{align*}
Define the deterministic variance matrix
\begin{align*}
V:=\sum_{i=1}^n\mathbb E[Y_i^2]\in\mathbb S^d.
\end{align*}
Each matrix $Y_i^2$ is positive semidefinite almost surely, so $V$ is positive semidefinite and $\sigma^2=\|V\|_{\mathrm{op}}=\lambda_{\max}(V)$.
If $L=0$, then $\|Y_i\|_{\mathrm{op}}=0$ almost surely for every $i$, hence $Y_i=0$ almost surely for every $i$. For $t=0$, the right-hand side is $2d$ by the zero-denominator convention in the theorem statement, so $\mathbb P(0\ge 0)=1\le 2d$. For $t>0$, we have $\mathbb P(0\ge t)=0$, so the desired inequality holds. If $\sigma^2=0$, then $V=0$. Since each $\mathbb E[Y_i^2]$ is positive semidefinite and
\begin{align*}
0=\operatorname{tr}V=\sum_{i=1}^n\mathbb E[\operatorname{tr}(Y_i^2)],
\end{align*}
we have $\operatorname{tr}(Y_i^2)=0$ almost surely for every $i$, hence $Y_i=0$ almost surely for every $i$. Thus the conclusion is again immediate, with the same convention handling $t=0$ and with probability $0$ when $t>0$. Hence it remains to prove the estimate when $L>0$ and $\sigma^2>0$.
[/step]
[step:Bound each matrix moment generating function by the variance contribution]
Define the scalar function $g:(0,3/L)\to(0,\infty)$ by
\begin{align*}
g(\theta)=\frac{\theta^2}{2(1-\theta L/3)}.
\end{align*}
Fix $\theta\in(0,3/L)$. For every $i\in\{1,\dots,n\}$,
\begin{align*}
\mathbb E[\exp(\theta Y_i)]\preceq \exp\left(g(\theta)\mathbb E[Y_i^2]\right),
\end{align*}
where $\preceq$ denotes the Loewner order on $\mathbb S^d$.
Indeed, for every real number $x$ with $|x|\le L$,
\begin{align*}
e^{\theta x}\le 1+\theta x+g(\theta)x^2.
\end{align*}
This is the scalar Bernstein exponential estimate; it follows from the power series for $e^{\theta x}$ and the bound
\begin{align*}
\frac{\theta^m x^m}{m!}\le \frac{\theta^m L^{m-2}x^2}{m!}
\end{align*}
for all $m\ge 2$ whenever $|x|\le L$, followed by
\begin{align*}
\sum_{m=2}^{\infty}\frac{\theta^mL^{m-2}}{m!}
\le \frac{\theta^2}{2}\sum_{m=0}^{\infty}\left(\frac{\theta L}{3}\right)^m
=\frac{\theta^2}{2(1-\theta L/3)}.
\end{align*}
Since $\|Y_i\|_{\mathrm{op}}\le L$ almost surely and $Y_i$ is symmetric almost surely, functional calculus applies this scalar inequality to every eigenvalue of $Y_i$ and gives
\begin{align*}
\exp(\theta Y_i)\preceq I+\theta Y_i+g(\theta)Y_i^2
\end{align*}
almost surely. Taking entrywise expectations preserves the Loewner order, and $\mathbb E[Y_i]=0$, so
\begin{align*}
\mathbb E[\exp(\theta Y_i)]
\preceq I+g(\theta)\mathbb E[Y_i^2].
\end{align*}
Since $g(\theta)\mathbb E[Y_i^2]$ is positive semidefinite, applying the scalar inequality $1+x\le e^x$ to its eigenvalues gives
\begin{align*}
I+g(\theta)\mathbb E[Y_i^2]\preceq \exp\left(g(\theta)\mathbb E[Y_i^2]\right).
\end{align*}
Combining the last two Loewner inequalities proves the claimed matrix moment generating function bound.
[guided]
The purpose of this step is to convert the almost sure operator norm bound into a usable exponential estimate. Define the scalar function $g:(0,3/L)\to(0,\infty)$ by
\begin{align*}
g(\theta)=\frac{\theta^2}{2(1-\theta L/3)}.
\end{align*}
Fix $\theta\in(0,3/L)$.
We first prove the scalar inequality that drives the whole argument. If $|x|\le L$, then the power series for the exponential gives
\begin{align*}
e^{\theta x}=1+\theta x+\sum_{m=2}^{\infty}\frac{\theta^m x^m}{m!}.
\end{align*}
For $m\ge 2$ and $|x|\le L$, we have $x^m\le |x|^m\le L^{m-2}x^2$, hence
\begin{align*}
\sum_{m=2}^{\infty}\frac{\theta^m x^m}{m!}
\le x^2\sum_{m=2}^{\infty}\frac{\theta^mL^{m-2}}{m!}.
\end{align*}
The Bernstein simplification replaces the exact exponential remainder by a geometric upper bound. Since $m!\ge 2\cdot 3^{m-2}$ for every $m\ge 2$, we get
\begin{align*}
\sum_{m=2}^{\infty}\frac{\theta^mL^{m-2}}{m!}
\le \frac{\theta^2}{2}\sum_{m=0}^{\infty}\left(\frac{\theta L}{3}\right)^m
=\frac{\theta^2}{2(1-\theta L/3)}
=g(\theta),
\end{align*}
where the geometric series converges because $\theta L/3<1$. Therefore
\begin{align*}
e^{\theta x}\le 1+\theta x+g(\theta)x^2
\end{align*}
for every $x\in[-L,L]$.
Now apply this scalar inequality to the eigenvalues of the symmetric random matrix $Y_i$. Since $\|Y_i\|_{\mathrm{op}}\le L$ almost surely, every eigenvalue of $Y_i$ belongs to $[-L,L]$ almost surely. Functional calculus therefore gives the Loewner-order inequality
\begin{align*}
\exp(\theta Y_i)\preceq I+\theta Y_i+g(\theta)Y_i^2
\end{align*}
almost surely. Taking expectations is legitimate because all entries are bounded random variables, and expectation preserves the Loewner order: if $A(\omega)\preceq B(\omega)$ almost surely, then $\mathbb E[A]\preceq\mathbb E[B]$. Since $Y_i$ is centered, $\mathbb E[Y_i]=0$, so
\begin{align*}
\mathbb E[\exp(\theta Y_i)]
\preceq I+g(\theta)\mathbb E[Y_i^2].
\end{align*}
Finally, $\mathbb E[Y_i^2]$ is positive semidefinite, because $Y_i^2$ is positive semidefinite almost surely. Hence $g(\theta)\mathbb E[Y_i^2]$ is positive semidefinite. Applying $1+x\le e^x$ to each eigenvalue gives
\begin{align*}
I+g(\theta)\mathbb E[Y_i^2]\preceq \exp\left(g(\theta)\mathbb E[Y_i^2]\right).
\end{align*}
Combining the two Loewner inequalities yields
\begin{align*}
\mathbb E[\exp(\theta Y_i)]\preceq \exp\left(g(\theta)\mathbb E[Y_i^2]\right).
\end{align*}
[/guided]
[/step]
[step:Combine the independent summands inside a trace exponential]
We use the [Lieb trace moment generating function estimate](https://doi.org/10.1007/s00440-011-0379-z), cited here as Tropp, *User-Friendly Tail Bounds for Sums of Random Matrices*, Theorem 3.6. We use the following Loewner-dominated corollary of that theorem: if $X_1,\dots,X_n$ are independent symmetric random matrices and $A_i$ are deterministic symmetric matrices such that
\begin{align*}
\mathbb E[\exp(X_i)]\preceq \exp(A_i)
\end{align*}
for every $i$, then
\begin{align*}
\mathbb E\left[\operatorname{tr}\exp\left(\sum_{i=1}^nX_i\right)\right]
\le \operatorname{tr}\exp\left(\sum_{i=1}^nA_i\right).
\end{align*}
Indeed, $\exp(X_i)$ is positive definite almost surely, so $\mathbb E[\exp(X_i)]$ is positive definite and its matrix logarithm is well-defined. Tropp's theorem is stated with $\log\mathbb E[\exp(X_i)]$ in place of $A_i$. From $\mathbb E[\exp(X_i)]\preceq\exp(A_i)$ and the operator monotonicity of the matrix logarithm on positive definite matrices, we get
\begin{align*}
\log\mathbb E[\exp(X_i)]\preceq A_i.
\end{align*}
The trace exponential map $B\mapsto\operatorname{tr}\exp(B)$ is monotone in the Loewner order on $\mathbb S^d$, so replacing $\log\mathbb E[\exp(X_i)]$ by the larger matrix $A_i$ preserves the upper bound. This is the finite-dimensional noncommutative trace moment generating function estimate obtained from Lieb's concavity theorem.
Apply this estimate with $X_i:=\theta Y_i$ and
\begin{align*}
A_i:=g(\theta)\mathbb E[Y_i^2]\in\mathbb S^d.
\end{align*}
The matrices $\theta Y_i$ are independent and symmetric almost surely, and the preceding step verifies the required moment generating function hypothesis. Therefore
\begin{align*}
\mathbb E[\operatorname{tr}\exp(\theta S)]
\le \operatorname{tr}\exp\left(g(\theta)V\right).
\end{align*}
Since $V$ is positive semidefinite and $\lambda_{\max}(V)=\sigma^2$, every eigenvalue of $g(\theta)V$ is at most $g(\theta)\sigma^2$. Thus
\begin{align*}
\operatorname{tr}\exp\left(g(\theta)V\right)
\le d\exp\left(g(\theta)\sigma^2\right).
\end{align*}
Consequently,
\begin{align*}
\mathbb E[\operatorname{tr}\exp(\theta S)]
\le d\exp\left(\frac{\theta^2\sigma^2}{2(1-\theta L/3)}\right).
\end{align*}
[guided]
This step is where independence is used in a genuinely matrix-valued way. We apply the [Lieb trace moment generating function estimate](https://doi.org/10.1007/s00440-011-0379-z) in the following Loewner-dominated form: if $X_1,\dots,X_n$ are independent symmetric random matrices and deterministic symmetric matrices $A_1,\dots,A_n$ satisfy
\begin{align*}
\mathbb E[\exp(X_i)]\preceq\exp(A_i)
\end{align*}
for every $i$, then
\begin{align*}
\mathbb E\left[\operatorname{tr}\exp\left(\sum_{i=1}^nX_i\right)\right]
\le \operatorname{tr}\exp\left(\sum_{i=1}^nA_i\right).
\end{align*}
The usual statement has $\log\mathbb E[\exp(X_i)]$ in place of $A_i$. This is compatible with the displayed form because $\exp(X_i)$ is positive definite almost surely, hence $\mathbb E[\exp(X_i)]$ is positive definite and its logarithm is defined. The matrix logarithm is operator monotone on positive definite matrices, and $B\mapsto\operatorname{tr}\exp(B)$ is monotone in the Loewner order on $\mathbb S^d$.
We apply the estimate with $X_i:=\theta Y_i$ and
\begin{align*}
A_i:=g(\theta)\mathbb E[Y_i^2]\in\mathbb S^d.
\end{align*}
The matrices $\theta Y_i$ are independent because the $Y_i$ are independent, and they are symmetric almost surely because each $Y_i$ is symmetric almost surely. The preceding step proves exactly the required hypothesis
\begin{align*}
\mathbb E[\exp(\theta Y_i)]\preceq\exp\left(g(\theta)\mathbb E[Y_i^2]\right).
\end{align*}
Therefore
\begin{align*}
\mathbb E[\operatorname{tr}\exp(\theta S)]
\le \operatorname{tr}\exp\left(g(\theta)V\right).
\end{align*}
Since $V$ is positive semidefinite and $\lambda_{\max}(V)=\sigma^2$, every eigenvalue of $g(\theta)V$ is at most $g(\theta)\sigma^2$. Summing the exponentials of the $d$ eigenvalues gives
\begin{align*}
\operatorname{tr}\exp\left(g(\theta)V\right)
\le d\exp\left(g(\theta)\sigma^2\right).
\end{align*}
Substituting the definition of $g(\theta)$ yields
\begin{align*}
\mathbb E[\operatorname{tr}\exp(\theta S)]
\le d\exp\left(\frac{\theta^2\sigma^2}{2(1-\theta L/3)}\right).
\end{align*}
[/guided]
[/step]
[step:Derive and optimise the upper tail bound for the largest eigenvalue]
For $t\ge 0$ and $\theta\in(0,3/L)$, let $\mathbb 1_{\{\lambda_{\max}(S)\ge t\}}:\Omega\to\{0,1\}$ denote the indicator function of the event $\{\omega\in\Omega:\lambda_{\max}(S(\omega))\ge t\}$. This event implies
\begin{align*}
e^{\theta\lambda_{\max}(S)}\ge e^{\theta t}.
\end{align*}
Since $\operatorname{tr}\exp(\theta S)=\sum_{j=1}^d e^{\theta\lambda_j(S)}\ge e^{\theta\lambda_{\max}(S)}$, where $\lambda_1(S),\dots,\lambda_d(S)$ are the eigenvalues of $S$, we obtain the pointwise inequality
\begin{align*}
\mathbb 1_{\{\lambda_{\max}(S)\ge t\}}
\le e^{-\theta t}\operatorname{tr}\exp(\theta S).
\end{align*}
Taking expectations and using the trace moment generating function bound gives
\begin{align*}
\mathbb P(\lambda_{\max}(S)\ge t)
\le d\exp\left(-\theta t+\frac{\theta^2\sigma^2}{2(1-\theta L/3)}\right).
\end{align*}
For $t>0$, choose
\begin{align*}
\theta:=\frac{t}{\sigma^2+Lt/3}.
\end{align*}
Because $L>0$ and $\sigma^2>0$, this choice satisfies $0<\theta<3/L$. With $a:=L/3$, we compute
\begin{align*}
1-a\theta
=1-\frac{at}{\sigma^2+at}
=\frac{\sigma^2}{\sigma^2+at}.
\end{align*}
Therefore
\begin{align*}
-\theta t+\frac{\theta^2\sigma^2}{2(1-a\theta)}
=-\frac{t^2}{\sigma^2+at}
+\frac{t^2}{(\sigma^2+at)^2}\cdot \frac{\sigma^2}{2}\cdot \frac{\sigma^2+at}{\sigma^2}.
\end{align*}
Simplifying the right-hand side gives
\begin{align*}
-\theta t+\frac{\theta^2\sigma^2}{2(1-a\theta)}
=-\frac{t^2}{2(\sigma^2+at)}.
\end{align*}
Since $a=L/3$, this is
\begin{align*}
-\theta t+\frac{\theta^2\sigma^2}{2(1-a\theta)}
=-\frac{t^2/2}{\sigma^2+Lt/3}.
\end{align*}
Hence, for every $t>0$,
\begin{align*}
\mathbb P(\lambda_{\max}(S)\ge t)
\le d\exp\left(-\frac{t^2/2}{\sigma^2+Lt/3}\right).
\end{align*}
For $t=0$, the same inequality holds because the left-hand side is at most $1$ and the right-hand side is $d\ge 1$.
[guided]
Fix $t\ge0$ and $\theta\in(0,3/L)$. Let $\mathbb 1_{\{\lambda_{\max}(S)\ge t\}}:\Omega\to\{0,1\}$ denote the indicator function of the event $\{\omega\in\Omega:\lambda_{\max}(S(\omega))\ge t\}$. This event forces
\begin{align*}
e^{\theta\lambda_{\max}(S)}\ge e^{\theta t},
\end{align*}
because the scalar exponential is increasing and $\theta>0$. If $\lambda_1(S),\dots,\lambda_d(S)$ are the eigenvalues of the symmetric matrix $S$, then
\begin{align*}
\operatorname{tr}\exp(\theta S)=\sum_{j=1}^d e^{\theta\lambda_j(S)}\ge e^{\theta\lambda_{\max}(S)}.
\end{align*}
Thus the indicator of the upper-tail event is pointwise bounded by
\begin{align*}
\mathbb 1_{\{\lambda_{\max}(S)\ge t\}}
\le e^{-\theta t}\operatorname{tr}\exp(\theta S).
\end{align*}
Taking expectations with respect to $\mathbb P$ and using the trace moment generating function bound from the previous step gives
\begin{align*}
\mathbb P(\lambda_{\max}(S)\ge t)
\le d\exp\left(-\theta t+\frac{\theta^2\sigma^2}{2(1-\theta L/3)}\right).
\end{align*}
For $t>0$, choose
\begin{align*}
\theta:=\frac{t}{\sigma^2+Lt/3}.
\end{align*}
Since $L>0$ and $\sigma^2>0$, this choice lies in $(0,3/L)$. Put $a:=L/3$. Then
\begin{align*}
1-a\theta
=1-\frac{at}{\sigma^2+at}
=\frac{\sigma^2}{\sigma^2+at}.
\end{align*}
Substitution gives
\begin{align*}
-\theta t+\frac{\theta^2\sigma^2}{2(1-a\theta)}
=-\frac{t^2}{\sigma^2+at}
+\frac{t^2}{(\sigma^2+at)^2}\cdot\frac{\sigma^2}{2}\cdot\frac{\sigma^2+at}{\sigma^2}
=-\frac{t^2}{2(\sigma^2+at)}.
\end{align*}
Since $a=L/3$, this becomes
\begin{align*}
\mathbb P(\lambda_{\max}(S)\ge t)
\le d\exp\left(-\frac{t^2/2}{\sigma^2+Lt/3}\right).
\end{align*}
For $t=0$, the same estimate holds because the left-hand side is at most $1$ and the right-hand side is $d\ge1$.
[/guided]
[/step]
[step:Apply the upper tail estimate to the negative sum and combine the tails]
Define the random matrix $\widetilde S:\Omega\to\mathbb S^d$ by
\begin{align*}
\widetilde S(\omega)=-S(\omega)=\sum_{i=1}^n(-Y_i(\omega)).
\end{align*}
The matrices $-Y_1,\dots,-Y_n$ are independent, centered, symmetric random matrices, and
\begin{align*}
\|-Y_i\|_{\mathrm{op}}=\|Y_i\|_{\mathrm{op}}\le L
\end{align*}
almost surely for every $i$. Their variance matrix is unchanged because $(-Y_i)^2=Y_i^2$, so the same value of $\sigma^2$ applies. Applying the largest-eigenvalue bound to $\widetilde S=-S$ gives
\begin{align*}
\mathbb P(\lambda_{\max}(-S)\ge t)
\le d\exp\left(-\frac{t^2/2}{\sigma^2+Lt/3}\right).
\end{align*}
For every symmetric matrix $A\in\mathbb S^d$,
\begin{align*}
\|A\|_{\mathrm{op}}=\max\{\lambda_{\max}(A),\lambda_{\max}(-A)\}.
\end{align*}
Thus
\begin{align*}
\left\{\|S\|_{\mathrm{op}}\ge t\right\}
\subset
\left\{\lambda_{\max}(S)\ge t\right\}
\cup
\left\{\lambda_{\max}(-S)\ge t\right\}.
\end{align*}
Using the union bound and the two tail estimates,
\begin{align*}
\mathbb P(\|S\|_{\mathrm{op}}\ge t)
\le \mathbb P(\lambda_{\max}(S)\ge t)+\mathbb P(\lambda_{\max}(-S)\ge t).
\end{align*}
The right-hand side is bounded by
\begin{align*}
2d\exp\left(-\frac{t^2/2}{\sigma^2+Lt/3}\right).
\end{align*}
Since $S=\sum_{i=1}^nY_i$, this is the asserted Matrix Bernstein inequality.
[guided]
The largest-eigenvalue estimate controls only the upper spectral tail. To control the operator norm of a symmetric matrix, we also need the lower spectral tail. Define the random matrix $\widetilde S:\Omega\to\mathbb S^d$ by
\begin{align*}
\widetilde S(\omega)=-S(\omega)=\sum_{i=1}^n(-Y_i(\omega)).
\end{align*}
The matrices $-Y_1,\dots,-Y_n$ remain independent because [measurable functions](/page/Measurable%20Functions) of independent random variables are independent. They are centered because $\mathbb E[-Y_i]=-\mathbb E[Y_i]=0$, and they are symmetric almost surely because each $Y_i$ is symmetric almost surely. Their operator norms satisfy
\begin{align*}
\|-Y_i\|_{\mathrm{op}}=\|Y_i\|_{\mathrm{op}}\le L
\end{align*}
almost surely, and their variance matrix is the same since $(-Y_i)^2=Y_i^2$. Therefore the upper-tail estimate applied to $\widetilde S=-S$ gives
\begin{align*}
\mathbb P(\lambda_{\max}(-S)\ge t)
\le d\exp\left(-\frac{t^2/2}{\sigma^2+Lt/3}\right).
\end{align*}
For every symmetric matrix $A\in\mathbb S^d$, the operator norm is the larger of the positive and negative spectral extremes:
\begin{align*}
\|A\|_{\mathrm{op}}=\max\{\lambda_{\max}(A),\lambda_{\max}(-A)\}.
\end{align*}
Hence
\begin{align*}
\left\{\|S\|_{\mathrm{op}}\ge t\right\}
\subset
\left\{\lambda_{\max}(S)\ge t\right\}
\cup
\left\{\lambda_{\max}(-S)\ge t\right\}.
\end{align*}
Applying the union bound gives
\begin{align*}
\mathbb P(\|S\|_{\mathrm{op}}\ge t)
\le \mathbb P(\lambda_{\max}(S)\ge t)+\mathbb P(\lambda_{\max}(-S)\ge t).
\end{align*}
Applying the two one-sided estimates to the right-hand side gives
\begin{align*}
\mathbb P(\|S\|_{\mathrm{op}}\ge t)
\le 2d\exp\left(-\frac{t^2/2}{\sigma^2+Lt/3}\right).
\end{align*}
Since $S=\sum_{i=1}^nY_i$, this is exactly the claimed operator-norm bound.
[/guided]
[/step]
