I will formulate my question first; below is the relevant background information and notation. I have asked this question on stack, but I think the odds are better that I actually get an answere here (but do correct me if I'm wrong to post it here).
WHY is the Shapovalov form on a Verma module symmetric?
The fact that it's orthogonal with respect to its weight decomposition simplifies things of course, so that we only have to calculate expressions of the form $$ \langle f^{n_1}_1\ldots f^{n_l}_lv_\lambda,f^{m_1}_1\ldots f^{m_l}_lv_\lambda\rangle $$ where $n,m\in\mathbb{Z}^l_{\ge0}$ are such that $\sum^l_{i=1}n_i\beta_i=\sum^l_{i=1}m_i\beta_i$. This seems rather hard to do because one has to commute elements of root spaces of non-simple roots. In principle, any element of the form $f^{n_1}_1\ldots f^{n_l}_l$ for $n\in\mathbb{Z}^l_{\ge0}$ can be expressed as a sum of elements of the form $f^{n_{i_1}}_{i_1}\ldots f^{n_{i_p}}_{i_p}$, where now $p\in\mathbb{Z}_{\ge1}$, and $i_1,\ldots,i_p\in\{1,\ldots,r\}$ (not necessarily distinct!) and $n_{i_1},\ldots,n_{i_l}\in\mathbb{Z}_{\ge1}$. This seems hopeless as well. In general, I don't feel like such a computation-based approach can be very helpful (let alone insightful). An elementary argument would be awesome.
Background
Let $L$ be a simple Lie algebra (over $\mathbb{C}$), let $L=L_-\oplus H\oplus L_+$ be a root space decomposition with $H$ a Cartan subalgebra of $L$ ($\dim(H)=r$ and $\dim(L_-)=\dim(L_+)=l$), let $\Delta=\{\alpha_1,\ldots,\alpha_r\}$ be a set of simple positive roots, and enumerate the positive roots as $\beta_1,\ldots,\beta_r$ such that the first $r$ are the simple onese. Denote by $B$ the Borel subalgebra $H\oplus L^+$. Let $\{e_1,\ldots,e_r,f_1,\ldots,f_r\}$ be a set of Chevalley generators, and let $\{e_{r+1},\ldots,e_l,f_{r+1},\ldots,f_l\}$ be elements of the root spaces corresponding to the positive, negative roots that are not simple.
Let $\lambda\in H^\ast$, denote by $\lambda^+\in B^\ast$ its extension to $B$ (by declaring it to be zero on $L_+$), and denote by $\mathbb{C}_\lambda$ the one-dimensional (over $\mathbb{C}$) $\mathfrak{U}(B)$-module defined by $xv_\lambda=\lambda^+(x)v$ for any $x\in B$ and for any $v\in\mathbb{C}_\lambda$, where $v$ is a nonzero element of $\mathbb{C}_\lambda$. The Verma module is then the induced module $M_\lambda:=\mathfrak{U}(L)\otimes_{\mathfrak{U}(B)}\mathbb{C}_\lambda$. Denote $1\otimes 1\in M_\lambda$ by $v_\lambda$.
The Verma module is a direct sum of weight spaces that are finite-dimensional, i.e. $$ M_\lambda=\bigoplus_{k\in\mathbb{Z}^r_{\ge0}}M^k_\lambda $$ where $M^k_\lambda=\{v\in M_\lambda\mid\forall\,h\in H:hv=(\lambda-\sum^r_{i=1}k_i\alpha_i)(h)v\}$ for any $k\in\mathbb{Z}^r_{\ge0}$. It is easy to see that, for any $k\in\mathbb{Z}^r_{\ge0}$, the weight space $M^k_\lambda$ is the linear span of $$ \{f^{n_1}_1\ldots f^{n_l}_l\mid n\in\mathbb{Z}^l_{\ge0}\,\wedge\,\sum^l_{i=1}n_i\beta_i=\sum^r_{i=1}k_i\alpha_i\}. $$ The restricted dual $M^\vee_\lambda$ of $M_\lambda$ corresponding to the Chevalley involution $\omega\colon L\to L$ defined by $$ \forall\,i\in\{1,\ldots,r\}\colon\omega(e_i)=-f_i\,\wedge\,\omega(f_i)=-e_i, $$ is as the $\mathfrak{U}(L)$-module whose underlying abelian group is $$ \bigoplus_{k\in\mathbb{Z}^r_{\ge0}}(M^k_\lambda)^\ast $$ and where the action of $\mathfrak{U}(L)$ is defined for each $k\in\mathbb{Z}^r_{\ge0}$ by setting $$ (xf)(v)=f(-\omega(x)v) $$ for all $x\in L$, for all $f\in(M^k_\lambda)^\ast$, and for each $v\in M^k_\lambda$. The unique vector $v^\ast_\lambda\in(M^0)^\ast_\lambda$ dual to $v_\lambda$ is an element of the weight space $(M^\vee_\lambda)^\lambda$ and it's annihilated by $L_+$, so we have a unique morphism $\phi\colon M_\lambda\to M^\vee_\lambda$ of left $\mathfrak{U}(L)$-modules such that $\phi(v_\lambda)=v^\ast_\lambda$. Now we define the Shapovalov form to be the bilinear form $$ \langle\cdot,\cdot\rangle:M_\lambda\times M_\lambda\to\mathbb{C},(v,w)\mapsto\phi(v)(w). $$ It has the property that $\langle M^k_\lambda,M^{k^\prime}_\lambda\rangle=\{0\}$ for all $k,k^\prime\in\mathbb{Z}^r_{\ge0}$ such that $k\neq k^\prime$. In addition, it satisfies $\langle v_\lambda,v_\lambda\rangle=1$ and $\langle xv,w\rangle=\langle v,-\omega(x)w\rangle$ for any $x\in L$ and for any $v,w\in M_\lambda$.