Let $A \in \mathbb{C}^{n \times n}$ and for a given $x \in \mathbb{C}^n$, the square gain of $A$ is $\lVert A x \rVert^2 / \lVert x \rVert^2$ where the norm is the usual Euclidean norm. I want to calculate the average gain over all possible $x \in \mathbb{C}^n$. Obviously, the gain is independent from the "length" of $x$, so without losing generality we can assume $\lVert x \rVert = 1$.
Now we need a method to parameterize all $\lVert x \rVert = 1$. The first thing that comes to my mind is using the singular value decomposition of $A$. Let $A = U \Sigma V^*$ be the SVD of $A$, so
$$x = \sum_{i=1}^n \alpha_i v_i ~~~~\text{such that}~~~~ \sum_{i=1}^n \alpha_i^2 = 1 $$
characterize all possible vectors, where $U = [u_1 ~~ \dots ~~ u_n], \Sigma = \operatorname{diag}\{\sigma_i\}$ and $V = [v_1 ~~ \dots ~~ v_n]$.
Let's first consider $n=2$. We can select $\alpha_1 := \cos \theta$ and $\alpha_2 := \sin \theta$. So, the average gain can be calculated as
$$\begin{align*} \frac{1}{2 \pi} \int_0^{2\pi} \lVert A v_1 \cos \theta + A v_2 \sin \theta \rVert^2 d\theta &= \frac{1}{2 \pi} \int_0^{2\pi} \lVert \sigma_1 u_1 \cos \theta + \sigma_2 u_2 \sin \theta \rVert^2 d\theta \\ &= \frac{1}{2 \pi} \int_0^{2\pi} \left( \sigma_1^2 \cos^2 \theta + \sigma_2^2 \sin^2 \theta \right) d\theta \\ &= \frac{\sigma_1^2 + \sigma_2^2}{2} \end{align*}$$
This result makes intuitive sense since the average square gain is the average of minimum and maximum square gains. So, we can hypotesize that the average square gain is $\lVert A \rVert_F^2 / n$ where $\lVert \cdot \rVert_F$ is the Frobenius norm. Now we can test this idea for $n=3$ using spherical coordinates. So let $\alpha_1 := \cos \theta, \alpha_2 := \sin \theta \cos \phi, \alpha_3 := \sin \theta \sin \phi$. Using the same steps we get
$$\begin{align*} &\frac{1}{2 \pi^2} \int_0^{\pi} \int_0^{2\pi} \lVert A v_1 \cos \theta + A v_2 \sin \theta \cos \phi + A v_3 \sin \theta \sin \phi \rVert^2 d\phi d\theta \\ &~~~~~~~~= \frac{1}{2 \pi^2} \int_0^{\pi} \int_0^{2\pi} \left( \sigma_1^2 \cos^2 \theta + \sigma_2^2 \sin^2 \theta \cos^2 \phi + \sigma_3^2 \sin^2 \theta \sin^2 \phi \right) d\phi d\theta \\ &~~~~~~~~= \frac{2\sigma_1^2 + \sigma_2^2 + \sigma_3^2}{4} \end{align*}$$
which is not what we expected. Coming back to $n=2$ case, we could have also selected $\alpha_1 = t$ and $\alpha_2 = \sqrt{1-t^2}$ where $t \in [0,1]$. In this case, we get
$$ \int_0^1 \left(\sigma_2^2 + (\sigma_1^2 - \sigma_2^2) t^2 \right) dt = \frac{\sigma_1^2 + 2 \sigma_2^2}{3} $$
which is different from our first attempt.
So, why are the results are different? Is there a "correct" way of selecting the parameters to get a "nice" formula like $\lVert A \rVert_F^2 / n$?
Equip $S^{n - 1} = \{x \in \mathbb{R}^n : |x| = 1\}$ with it's usual surface measure, normalized to be a probability measure. Then you seek $\int_{S^{n - 1}}|Ax|^2\,dS(x)$.
As for selecting parameters, any parameterization of $S^{n - 1}$ will work. You have to apply the formula for the surface integral though:
If $\phi : O \to U \subset {S}^{n - 1}$ is a $C^1$ parameterization and $f : S^{n - 1} \to \mathbb{R}$ vanishes off $U$, then $$\int_{U}f(x)\,dS(x) = \frac{1}{\text{Area}(S^{n - 1})}\int_{O}f(\phi(x))\sqrt{\det D\phi(x)^T D\phi(x)}\,dx.$$ I think you forgot the $\sqrt{ \det D\phi(x)^T D\phi(x)}$. Luckily, this problem is simple enough that we don't need to use coordinates.
Write the SVD of $A$ as $A = UDV^T$. Then $|Ax|^2 = |DV^Tx|^2$. So by invariance of the measure under orthogonal transformations, $$\int_{S^{n - 1}}|Ax|^2\,dS(x) = \int_{S^{n - 1}}|DV^Tx|^2\,dS(x) = \int_{S^{n - 1}}|Dx|^2\,dS(x).$$ Now note $|Dx|^2 = (D^2x, x) = \sigma_1^2|x_1|^2 + \dots + \sigma_n^2|x_n|^2$. Hence $$\int_{S^{n - 1}}|Dx|^2\,dS(x) = \sum_{i = 1}^{n}\sigma_i^2\int_{S^{n - 1}}|x_i|^2\,dS(x).$$ Note that $a_i = \int_{S^{n - 1}}|x_i|^2\,dS(x)$ is independent of $i$ by invariance of the measure under orthogonal transformations. Moreover, $$a_1 + \dots + a_n = \int_{S^{n - 1}}|x|^2\,dS(x) = 1.$$ Hence $a_i = \frac{1}{n}$. In conclusion, $$\int_{S^{n - 1}}|Ax|^2\,dS(x) = \frac{\sigma_1^2 + \dots + \sigma_n^2}{n}.$$
Edit: Here is why $a_i = \int_{S^{n - 1}}|x_i|^2 dS(x)$ is independent of $i$. Fix $i$. Note that \begin{align} \int_{S^{n - 1}}|x_i|^2\,dS(x) &= \int_{S^{n - 1}}|(x, e_i)|^2\,dS(x). \end{align} Now pick an orthogonal transformation $R$ such that $Re_i = e_1$. Since $(x, e_i) = (Rx, Re_i) = (Rx, e_1)$, it follows that $$\int_{S^{n - 1}}|x_i|^2\,dS(x) = \int_{S^{n - 1}}|(Rx, e_1)|^2\,dS(x) = \int_{S^{n - 1}}|(x, e_1)|^2\,dS(x),$$ Where the last equality follows from invariance of the measure under orthogonal transformations. So $a_i = a_1$.
The definition of the "surface measure" I am using here is the one that works for any $C^1$ $m$-dimensional surface $M \subset \mathbb{R}^k$. The measure $\mu$ is a Borel measure on $M$ defined in such a way that if $f : M \to [0, \infty]$ is a measurable function vanishing off the image $U \subset M$ of a coordinate parameterization $\phi : O \to U$, $O$ open in $\mathbb{R}^m$, then $$\int_{M}f\,d\mu = \int_{U}f\,d\mu = \int_{O}f(\phi(x))\sqrt{g(x)}\,dx, \hspace{20pt} g(x) = \det D\phi(x)^T D\phi(x).$$ Here I used it for $M = S^{n - 1}$ (I also normalized it to $1$). From this you can prove the assertion that the measure on $S^{n - 1}$ is invariant under orthogonal transformations, that is, if $R^TR = I$, and $f : S^{n - 1} \to [0, \infty]$ is measurable, then $$\int_{S^{n - 1}}f(Rx)\,d\mu(x) = \int_{S^{n - 1}}f(x)\,d\mu(x).$$ A way to prove this is to first prove it for $f$ vanishing off of a coordinate patch, in which case it is just a computation. Then the fact that $M$ is always a countable union of such patches implies the general case.
The above measure is also defined on a general Riemannian manifold $M$, in which case $D\phi(x)^T D\phi(x)$ is replaced with the coordinate matrix $G(x)$ of the metric tensor of $M$ (the above definition is the special case when the metric tensor is the dot product).