In McDuff and Salamon's Intro. to Symplectic Topology (2nd edition), there's a proof that $U(n)$ is a maximal compact subgroup of $Sp(2n)$ which I'm trying to understand. The proof uses the Haar measure to average over matrices in a way that's unclear to me. The statement I'm having difficulties with is the following: Let $G\subset Sp(2n)$ be a compact subgroup. There exist a symmetric and positive definite matrix $P\in Sp(2n)$ such that $$\Psi^TP\Psi=P, \forall\Psi\in G$$ According to the authors, such a matrix can be obtained by averaging the matrices $\Psi^T\Psi$ over $\Psi\in G$ using the Haar measure $C(G,\mathbb R)\rightarrow \mathbb R$ for a compact Lie group.
So, the question is, in its broadest sense, how to understand the averaging process here? (Bear with me, clearer question follows)
Here are some details on my progress in making sense of this before submitting this question: So far I have came across these Notes on Compact Lie Groups by Salamon which detail the construction of said Haar measure. This does not explain how to obtain a measure, in the measure theoretic sense, from the functional built in Salamon's notes. So, further searching lead me to Royden's Real Analysis (3rd edition), where he explains in detail how to build a measure out of a functional (Daniell Integral) in such a way that the integral over the measure agrees with the functional. All of this is fine and well, but ultimately remains in the realm of functionals, where I'm looking for a process that'll result in the matrix $P$ as above.
In a more intuitive manner, what I gather should happen is somewhere along these lines: Let $M:C(G,\mathbb R)\rightarrow\mathbb R$ be a Haar measure following Salamon's notations, and let $\mu$ be the measure corresponding to the functional $M$, constructed like in Royden. Put $P=\int_G \Phi^T\Phi\operatorname{d}\mu$, then we can calculate that $$\Psi^TP\Psi=\Psi^T\int_G\Phi^T\Phi\operatorname{d}\mu\Psi=\int_G\Psi^T\Phi^T\Phi\Psi\operatorname{d}\mu=\int_G\left(\Phi\Psi\right)^T\Phi\Psi\operatorname{d}\mu=\int_G \Phi^T\Phi\operatorname{d}\mu=P$$where the second to last transition follows from $\mu$ being right invariant.
So, now my question boils down to this - In light of all that's said, seeing as the Daniell integral $M$ assigns a real number to a function, where we need to assign a matrix to a function, how can $\int_G\Phi^T\Phi\operatorname{d}\mu$ be understood?
Also, I'm guessing the second transition in the calculation follows from some means of limiting and applying the equation to the limits. If I'm mistaken, I would very much appreciate clarification on that as well, assuming the entire maneuver makes sense.
Note: I am aware of an errata for this edition which fixes the error here, namely, $P$ is not guaranteed to be symplectic, but that's besides the point.