What is the Hessian of the spectral norm?

Let us arrange all eigenvalues (and corresponding eigenvectors) in ascending order: $|\lambda_1| \leq \cdots |\lambda_{n-1}| < |\lambda_n|$. Then the elements of the desired Hessian are: $$\frac{\partial^2 \|A\|_2}{\partial A_{kl}\partial A_{ij}} = \frac{\partial^2 \lambda_{n}}{\partial A_{kl}\partial A_{ij}} = [\frac{\partial u^T_{n}}{\partial A_{kl}}]_i [u_{n}]_j + [u^T_{n}]_i [\frac{\partial u_{n}}{\partial A_{kl}}]_j$$

Here we use the fact that $A$ is diagonalizable, i.e., the set of eigenvectors $u_i, \ i=1,\cdots, n$ form a orthogonal basis. The basis is used to decompose $\frac{\partial u_n}{\partial A_{kl}} = \sum_{m=1}^n c_m u_m$. To find the coefficients $c_m \in \mathbb{R}$ one could follow the perturbation analysis (https://en.wikipedia.org/wiki/Eigenvalue_perturbation).

We just give the results: $$\frac{\partial u_{n}}{\partial A_{kl}} = \sum_{m=1}^{n-1} \frac{[u_m]_k[u_n]_l}{\lambda_n - \lambda_m} u_m$$ $$\frac{\partial u^T_{n}}{\partial A_{kl}} = \sum_{m=1}^{n-1} \frac{[u_m]_l[u_n]_k}{\lambda_n - \lambda_m} u_m$$ And the Hessian elements are: $$\frac{\partial^2 \|A\|_2}{\partial A_{kl}\partial A_{ij}} = \sum_{m=1}^{n-1} \frac{[u_n]_k[u_n]_j[u_m]_l[u_m]_i + [u_n]_l[u_n]_i[u_m]_k[u_m]_j}{ \lambda_n - \lambda_m} $$