Let $f:\mathrm{Mat}_n(\mathbb{C})\to\mathbb{C}$ be some function and let us suppose we want to make a change of variables in the integral $$ \int_{A\in \mathrm{Mat}_n(\mathbb{C})}f(A)\mathrm{d}{A} $$ from $A$ to $|A| U$, i.e., the polar decomposition of $A$, where $|A|\equiv\sqrt{A^\ast A}$ and $U$ is the unique partial isometry with kernel equal to that of $A$ (there is a theorem saying it exists).
What is the Jacobian matrix of the transformation $A \mapsto (|A|, U$)? I.e., what is $J$ such that the following equation holds: $$ \int_{A\in \mathrm{Mat}_n(\mathbb{C})}f(A)\mathrm{d}{A} = \int_{P\geq0,U^\ast U\,\mathrm{idempotent}}f(P U)|\det(J(P,U))|\mathrm{d}{P}\mathrm{d}{U}$$
I tried to calcualte it but I'm not getting anything simple. In particular, I've written $A = A_R + i A_I$ with $A_R = \frac{1}{2}(A+A^\ast); A_I = \frac{1}{2i}(A-A^\ast)$ so that $A$ is parametrized by two self-adjoint matrices. In turn, we may write $|A| = \exp(H_1) ; U = \exp(i H_2)$ for two self-adjoint matrices $H_1,H_2$ (assuming for a moment that $A$ is invertible so that $U$ is actually unitary). Hence we want to calculate the Jacobian of the transformation $(H_1,H_2)\mapsto(A_R,A_I)$ from $\mathrm{Herm}_n(\mathbb{C})^2\to \mathrm{Herm}_n(\mathbb{C})^2$.
This, however, starts to get ugly, with the differential of the exponential map for example being given by functional calculus of the adjoint super operator (https://en.wikipedia.org/wiki/Derivative_of_the_exponential_map) and having to use the determinant of a block matrix formula.
Is there an easier way out?
Possible solution:
In Edelman's PhD thesis there are given Jacobians to get from a matrix A to its LQ decomposition, and from its LQ decomposition to its Cholesky decomposition (Theorem 3.1). This possibly solves the problem as follows:
\begin{align} \int_{A\in \mathrm{Mat}_n(\mathbb{C})}f(A)\mathrm{d}{A} &= \int_{L\text{ lower triangular},\,U\text{ unitary}} f(LU)\prod_{i=1}^{n}L_{ii}^{2n-2i+1}\mathrm{d}{L}\mathrm{d}{U} \\ &=2^{-n}\int_{P\geq0,\,U\text{ unitary}} f(\sqrt{P}U)\mathrm{d}{P}\mathrm{d}{U}\\&=2^{-n}\int_{P\geq0,\,U\text{ unitary}} f(PU)|\det(P\otimes I+I\otimes P^\ast)|^2\mathrm{d}{P}\mathrm{d}{U}\\&=2^{-n}\int_{P\geq0,\,U\text{ unitary}} f(PU)\prod_{1\leq i,j\leq n}(\lambda_i(P)+\lambda_j(P))^2\mathrm{d}{P}\mathrm{d}{U}\end{align}
with the usual abuse of notation that $\mathrm{d}{L}$ integrates only over the non-zero elements of $L$, $\mathrm{d}{U}$ is the volume element within the unitary group, and $\mathrm{d}{P}$ the volume element on self-adjoint matrices (so only $n$ real and $\frac{1}{2}n(n-1)$ complex matrix elements). $\lambda_j(P)$ is the $j$th eigenvalue of the matrix $P$.
Remaining question: Why is the LQ-decomposition change of variables valid for complex matrices? A complex unitary $n\times n$ matrix is $n^2$ real parameters, whereas a lower triangular matrix is $n(n+1)$ real parameters. On the other hand, a complex matrix is $2n^2$ real parameters, so there seem to be $n$ real parameters too many in this decomposition? (This is not a problem if matrices have real entries). Note that for the Cholesky decomposition this is not an issue since then the lower triangular matrix has positive entries on its diagonal.
Could it be possible to make an LQ decomposition for complex matrices where the lower triangular has positive entries on the diagonal? Is this what Edelman is referring to?
Unfortunately, precisely for the complex LQ decomposition he does not give a reference nor a proof.
First some generalities:
For this statement check any book that contains a treatment of the Haar measure. For example it is Theorem 2.49 in Folland - A Course in Abstract Harmonic Analysis.
Now we apply this to your situation.
The first comment is that the non-invertible matrices are a Lebesgue nullset in $M_{n\times n}(\Bbb C)$, hence we may make the domain of integration smaller so that we are integrating over $GL_n(\Bbb C)$. Now $GL_n(\Bbb C)$ is a group and as such there is a Haar measure on it, to be precise we choose the normalisation so that the measure is $$dg= \frac{\prod_{ij}dg_{ij}}{|\det(g)|^{2n}},$$ hence $$\int_{M_{n\times n}}f(A)dA = \int_{GL_n} f(g) |\det(g)|^{2n}dg.$$ Next note that the unitaries are a compact subgroup of $GL_n$ and that $GL_n$ is unimodular, meaning $\Delta G=1$. Further each class $GL_n/U(n)$ has a unique representant being a positive matrix, this a restatement of the theorem you cite in your question. As such by the above discussion you have $$\int_{M_{n\times n}}f(A)\,dA=\int_{GL_n}f(g)|\det(g)|^{2n}\,dg = \int_{GL_n/U(n)} d[p]\int_{U(n)}du\, f(p\cdot u)\det(p)^{2n}$$ Now we are almost finished. What is left to do is to relate the two measures on the right-hand side with the measures we are interested in. The easier of the two is the integral over $U(n)$.
It is not entirely clear to me what measure you are using on $U(n)$ in your question, but I see only two possible definitions and they are both the same (up to a constant). On the one hand you have the Haar measure on $U(n)$, on the other hand $M_{n\times n}$ is an euclidean vector space and the Riemannian metric then restricts to a metric on the sub-manifold $U(n)$, which will give you a volume form on $U(n)$. However the scalar product on $M_{n\times n}$ is given by $\langle A, B\rangle = \mathrm{Tr}(A^* B)$, so multiplication with $U(n)$ preserves this scalar product and $U(n)$ acts by isometries on $M_{n\times n}$. In particular $U(n)$ acts by isometries on the induced metric on $U(n)$ and as such preservers the volume form. This means it must be equal, up to a constant, to the Haar measure on $U(n)$. Googling reveals the value of the constant.
Now what about $GL_n/U(n)$? As noted this can be identified with the set of strictly positive matrices, which are an open cone in $\mathrm{Herm}_{n\times n}$ and as such you may use the Lebesgue measure of that vector space to integrate things. We should compare this measure to the one our theorem above gives.
The good thing is that our theorem gives uniqueness of the measure on $GL_n/U(n)$. So we just need to find a measure on the positive matrices which is invariant under the (correct) action of $GL_n$ and then tune the constants. I haven't done this calculation, but my guess is that what you get must be $$d[p]= \frac{d\lambda}{\det(p)^{n}}$$
Where $d\lambda$ is the Lebesgue measure on hermitian $n\times n$ matrices. This would give you an end-result: $$\int_{M_{n\times n}}f(A)dA = \int_{\mathrm{Pos}_{n\times n}}dP\int_{U(n)} dU f(PU)\det(p)^n \cdot \mathrm{Const.}$$ The constant is $\frac1{\mathrm{Vol}(U(n))}$ times the normalisation constant of $G/K$. You can test this formula and extract the constants for example by integrating things over the unit ball.