Let $|\psi_\rho\rangle\in\mathcal S(\mathcal X\otimes\mathcal Y)$ be a purification of $\rho\in\mathrm D(\mathcal X)$. This means that
$$\operatorname{Tr}_{\mathcal Y}|\psi_\rho\rangle\!\langle\psi_\rho|=\rho
\Longleftrightarrow \psi_\rho\psi_\rho^\dagger=\rho
\Longleftrightarrow \psi_\rho=\sqrt\rho U,$$
where $\psi_\rho\in\mathrm{Lin}(\mathcal Y,\mathcal X)$ is the linear operator with matrix elements $(\psi_\rho)_{ij}=\langle i,j|\psi_\rho\rangle$, and in the last step we took the polar decomposition of $\psi_\rho$, so that $U\in\mathrm U(\mathcal Y,\mathcal X)$ is any partial isometry such that ${\rm range}(U)={\rm supp}(\rho)$.
Note that this also means that $\sqrt\rho UU^\dagger\sqrt\rho=\rho$, i.e. $U^\dagger$ acts as an isometry at least on the support of $\rho$. In fact, we could simply assume $UU^\dagger=I$ without losing anything, because the action of $U^\dagger$ outside of the support of $\rho$ is never relevant to this discussion. Hence why the result can be perhaps more simply stated saying that $\psi_\rho=\sqrt\rho V^\dagger$ for some isometry $V$.
You can equivalently write the above saying that the purifications of $\rho$ have the form
$$|\psi_\rho\rangle = \operatorname{vec}(\sqrt\rho U)$$
for any $\mathcal Y$ such that $\dim(\mathcal Y)\ge\operatorname{rank}(\rho)$, and any partial isometry $U\in\mathrm U(\mathcal Y,\mathcal X)$ with range equal to that of $\rho$.
I personally think the easier way to put this is to say that, if $\rho$ has eigendecomposition $\rho=\sum_i p_i|u_i\rangle\!\langle u_i|$, then its purifications are all and only the vectors of the form
$$|\psi_\rho\rangle = \sum_i \sqrt{p_i} |u_i, v_i\rangle$$
for any orthonormal set of vectors $\{|v_j\rangle\}_j$.
Note that for this to be possible, you need $\operatorname{rank}(\psi_\rho)=\operatorname{rank}(\rho)$, which is only possible if $\dim(\mathcal Y)\ge\operatorname{rank}(\rho)$. This is the only assumption that needs to be made on $\mathcal Y$.
Now, back to Uhlmann. Clearly, to take an inner product between purifications, you need them to be defined in spaces of the same dimensions. All finite-dimensional spaces are isomorphic, so it's not really meaningful to ask anything about the structure of such spaces aside from their dimensions. We can thus take some $|\psi_\rho\rangle$ and $|\psi_\sigma\rangle$ defined over spaces of the same dimensions, and then notice that
$$\langle \psi_\rho,\psi_\sigma\rangle = \operatorname{Tr}(VU^\dagger \sqrt\rho\sqrt\sigma),$$
for some $U,V\in\mathrm U(\mathcal Y,\mathcal X)$, for some $\mathcal Y$ with $\dim(\mathcal Y)\ge \min(\operatorname{rank}(\rho),\operatorname{rank}(\sigma))$.
The conclusion follows from here using standard matrix inequalities, discussed e.g. here.
More specifically, we get
$$\lvert \langle \psi_\rho,\psi_\sigma\rangle\rvert
\le \|VU^\dagger\|_{\rm op}{\rm Tr}|\sqrt\rho\sqrt\sigma|
\le {\rm Tr}|\sqrt\rho\sqrt\sigma|,$$
where the maximum is achieved when $VU^\dagger$ are unitary. Note that here $U,V$ are defined as surjective partial isometries, and thus $U^\dagger$ is an isometry. Their product is therefore also an isometry (and also a unitary, if the purification spaces of $\rho$ and $\sigma$ coincide, which we can however always assume without loss of generality) iff $U$ and $V$ have identical support.
Toy example
Purifications of $\rho$ —
Let's consider a concrete toy example.
Say $\rho=I/2$ be a maximally mixed qubit state.
Its purifications are then states of the form $\Psi=\frac12 V^\dagger$ for isometries $V$. Examples of such purifications might be
$$|\Psi_1\rangle= \frac1{\sqrt2} (|00\rangle+|11\rangle),
\quad\text{using}\quad V_1=I, \\
|\Psi_2\rangle = \frac1{\sqrt2}(|01\rangle+|10\rangle),
\quad\text{using}\quad V_2=X, \\
|\Psi_3\rangle = \frac{|0\rangle(|0\rangle+|1\rangle+|2\rangle)
+|1\rangle(|0\rangle+\omega_3|1\rangle+\omega_3^2|2\rangle}{\sqrt6},
\, V_3=\frac1{\sqrt3}\begin{pmatrix}1&1\\1&\omega_3\\1&\omega_3^2\end{pmatrix},
$$
where $\omega_3\equiv\exp(2\pi i/3)$.
Note that the support of all these $V$ is the same as that of $\rho$. And you can always assume this without loss of generality.
Purifications of $\sigma$ —
Let's say $\sigma=\operatorname{diag}(1/3,2/3)$. Using the same isometies as above would then give the purifications
$$|\Phi_1\rangle= \frac1{\sqrt3} (|00\rangle+\sqrt2|11\rangle),
\quad\text{using}\quad V_1=I, \\
|\Phi_2\rangle = \frac1{\sqrt3}(|01\rangle+\sqrt2|10\rangle),
\quad\text{using}\quad V_2=X, \\
|\Phi_3\rangle = \frac{|0\rangle(|0\rangle+|1\rangle+|2\rangle)
+\sqrt2|1\rangle(|0\rangle+\omega_3|1\rangle+\omega_3^2|2\rangle}{3},
\, V_3=\frac1{\sqrt3}\begin{pmatrix}1&1\\1&\omega_3\\1&\omega_3^2\end{pmatrix}.
$$
Inner products of purifications —
Let's try something like
$$\langle\Psi_1|\Phi_2\rangle
= \operatorname{tr}[(\sqrt\rho V_1^\dagger)^\dagger\sqrt\sigma V_2^\dagger]
= \frac{1}{\sqrt2}\operatorname{tr}[\sqrt\sigma X]=0.$$
Well, that obviously isn't what we expect, because we know the fidelity between $\rho$ and $\sigma$ isn't zero. If we instead try using $V_1$ for both states, we get instead
$$\langle\Psi_1|\Phi_1\rangle
= \operatorname{tr}[\sqrt\rho\sqrt\sigma]
= \frac{1+\sqrt2}{\sqrt6}.$$
And this is the maximum because it's trivial to check that in this case $\sqrt\rho\sqrt\sigma=|\sqrt\rho\sqrt\sigma|$.
Note that we could have equivalently used the other purifications, as $\langle\Psi_1|\Phi_1\rangle=\langle\Psi_2|\Phi_2\rangle=\langle\Psi_3|\Phi_3\rangle$. This happens because, as discussed before, the only thing that matters to the overlaps is the product $V_1^\dagger V_2$ for the two isometries.
Slightly less toyish examples
Example with less trivial states —
The example above is kinda trivial because $\rho$ and $\sigma$ commute (and in fact, $\rho$ commutes with everything). For a less trivial analysis consider
$$
\rho' \equiv p \mathbb{P}_0 + (1-p) \mathbb{P}_1,\,\, p\neq1/2,\qquad
\sigma' \equiv q \mathbb{P}_+ + (1-q)\mathbb{P}_-,\,\, q\neq1/2,
$$
using the notation $\mathbb{P}_\psi\equiv |\psi\rangle\!\langle\psi|$.
These ensure $[\rho',\sigma']\neq0$.
Defining purifications in analogy of the calculations above in terms of the isometries $V_i$, we have for example
$$\langle\Psi_1|\Phi_1\rangle
= \operatorname{tr}[\sqrt{\rho'}\sqrt{\sigma'}]= \frac{(\sqrt p+\sqrt{1-p})(\sqrt q+\sqrt{1-q})}{2},$$
but this time this isn't the max overlap.
After some calculation, you indeed find that
$$\operatorname{tr}|\sqrt{\rho'}\sqrt{\sigma'}|
=\frac{\sqrt{1+\sqrt{u}}+\sqrt{1-\sqrt{u}}}{2},
\quad u\equiv 1-16 pq(1-p)(1-q),$$
and you can verify numerically that this is always larger than the other result. I couldn't find any simple expression for the purifications/isometries that produce this fidelity.
Pure states —
If $\rho=\mathbb{P}_\psi$ and $\sigma=\mathbb{P}_\phi$ are pure states, we of course know that the fidelity is just $|\langle\psi|\phi\rangle|$. It might be tempting to think that we get this following the argument above by simply having the unitary as $U=I$, and thus $F=\operatorname{tr}[\sqrt\rho\sqrt\sigma]$. This is however not so:
$$\operatorname{tr}[\sqrt\rho\sqrt\sigma]=|\langle\psi|\phi\rangle|^2,$$
and we should note that the square here is a problem, because the formula $\operatorname{tr}|\sqrt\rho\sqrt\sigma|$ is for the "square root" fidelity. You can also see that this is wrong because the result is trivially reachable reasoning in terms of purifications: the "purifications" of $\mathbb{P}_\psi$ and $\mathbb{P}_\phi$ are just $|\psi\rangle$ and $|\phi\rangle$, and thus their overlap is just $|\langle\psi|\phi\rangle|$. The apparent discrepancy happens because for pure states $|\sqrt\rho\sqrt\sigma|\neq \sqrt\rho\sqrt\sigma$.
Instead, one has
$$|\sqrt\rho\sqrt\sigma|=|\mathbb{P}_\psi\mathbb{P}_\phi|
= |\langle\psi|\phi\rangle| \sqrt{\mathbb{P}_\psi}
= |\langle\psi|\phi\rangle| \mathbb{P}_\psi,$$
whose trace is then clearly just $|\langle\psi|\phi\rangle|$, as expected.
