[proofplan]
We use the spectral theorem for normal matrices, which guarantees the existence of a unitary diagonalisation $A = UDU^*$ with $D = \operatorname{diag}(\lambda_1, \ldots, \lambda_M)$. Expanding an arbitrary vector $u$ in the orthonormal eigenbasis and computing $\|Au\|_2^2$ using orthonormality reduces the operator norm to $\max_i |\lambda_i|$. The upper bound $\|A\|_2 \leq \rho(A)$ follows from bounding each $|\lambda_i|$ by $\rho(A)$, and the matching lower bound $\|A\|_2 \geq \rho(A)$ is attained by choosing $u$ to be the eigenvector corresponding to the eigenvalue of largest modulus.
[/proofplan]
[step:Diagonalise $A$ via the spectral theorem for normal matrices]
Since $A \in \mathbb{C}^{M \times M}$ is normal ($A^*A = AA^*$), the spectral theorem for normal matrices guarantees the existence of a unitary matrix $U \in \mathbb{C}^{M \times M}$ (so $U^*U = UU^* = I$) and a diagonal matrix $D = \operatorname{diag}(\lambda_1, \ldots, \lambda_M)$ such that
\begin{align*}
A = UDU^*.
\end{align*}
The columns $q_1, \ldots, q_M$ of $U$ form an orthonormal eigenbasis for $\mathbb{C}^M$: $Aq_i = \lambda_i q_i$ and $(q_i, q_j)_{\ell^2} = \delta_{ij}$.
[/step]
[step:Expand $\|Au\|_2$ in the eigenbasis and reduce to $\max_i |\lambda_i|$]
Let $u \in \mathbb{C}^M \setminus \{0\}$. Since $(q_1, \ldots, q_M)$ is an orthonormal basis, we write
\begin{align*}
u = \sum_{i=1}^{M} a_i q_i, \qquad a_i = (u, q_i)_{\ell^2}.
\end{align*}
Applying $A$:
\begin{align*}
Au = \sum_{i=1}^{M} a_i Aq_i = \sum_{i=1}^{M} \lambda_i a_i q_i.
\end{align*}
By the orthonormality of $(q_i)$, the Pythagorean theorem gives
\begin{align*}
\|Au\|_2^2 = \sum_{i=1}^{M} |\lambda_i|^2 |a_i|^2, \qquad \|u\|_2^2 = \sum_{i=1}^{M} |a_i|^2.
\end{align*}
Therefore
\begin{align*}
\frac{\|Au\|_2^2}{\|u\|_2^2} = \frac{\sum_{i=1}^{M} |\lambda_i|^2 |a_i|^2}{\sum_{i=1}^{M} |a_i|^2}.
\end{align*}
This is a weighted average of $|\lambda_1|^2, \ldots, |\lambda_M|^2$ with non-negative weights $|a_i|^2 / \sum_j |a_j|^2$ summing to $1$. A weighted average of non-negative reals is bounded above by the maximum and below by the minimum:
\begin{align*}
\min_i |\lambda_i|^2 \leq \frac{\|Au\|_2^2}{\|u\|_2^2} \leq \max_i |\lambda_i|^2 = \rho(A)^2.
\end{align*}
Taking the supremum over all $u \neq 0$ on the left-hand side:
\begin{align*}
\|A\|_2 = \sup_{u \neq 0} \frac{\|Au\|_2}{\|u\|_2} \leq \rho(A).
\end{align*}
[guided]
The key computation uses the orthonormality of the eigenbasis. Once we expand $u = \sum_i a_i q_i$ and apply $A$, the cross terms vanish:
\begin{align*}
\|Au\|_2^2 = \left(\sum_i \lambda_i a_i q_i, \sum_j \lambda_j a_j q_j\right)_{\ell^2} = \sum_{i,j} \lambda_i \overline{\lambda_j} a_i \overline{a_j} (q_i, q_j)_{\ell^2} = \sum_i |\lambda_i|^2 |a_i|^2.
\end{align*}
The cross terms with $i \neq j$ vanish because $(q_i, q_j)_{\ell^2} = 0$. This is precisely the Pythagorean theorem for orthonormal expansions — it converts the matrix norm computation into a scalar problem.
The ratio $\|Au\|_2^2 / \|u\|_2^2 = \sum_i |\lambda_i|^2 |a_i|^2 / \sum_i |a_i|^2$ is a convex combination of the values $|\lambda_i|^2$. A convex combination is bounded by its largest term, giving the upper bound $\leq \rho(A)^2$.
Why does normality matter? For a non-normal matrix, the eigenvectors (if they exist) need not be orthogonal. Without orthogonality, the cross terms in $\|Au\|_2^2$ do not vanish, and the ratio $\|Au\|_2 / \|u\|_2$ can exceed $\rho(A)$. This is the phenomenon of non-normal amplification.
[/guided]
[/step]
[step:Attain equality by choosing $u$ as the eigenvector with largest eigenvalue modulus]
Let $i_0 \in \{1, \ldots, M\}$ be an index such that $|\lambda_{i_0}| = \rho(A)$. Choose $u = q_{i_0}$. Then $a_{i_0} = 1$ and $a_j = 0$ for $j \neq i_0$, so
\begin{align*}
\frac{\|Aq_{i_0}\|_2}{\|q_{i_0}\|_2} = \frac{|\lambda_{i_0}| \cdot 1}{1} = |\lambda_{i_0}| = \rho(A).
\end{align*}
Therefore $\|A\|_2 \geq \rho(A)$. Combining with the upper bound from the previous step:
\begin{align*}
\|A\|_2 = \rho(A).
\end{align*}
[/step]