Yes. Christian Remling's comment in essence already gave an affirmative answer, and this question is a bit old, but I think it's worthwhile to write a proof that does not require reducing to the separable situation.

I encountered the same question myself and have not found any reference---but with the polar decomposition and spectral theorems, the result follows with little fuss; this leads me to be believe it is almost certainly known, if not written down somewhere on some ancient scroll of functional analysis. (But the fact that I haven't seen it written down in general form compels me to put it out into the internet ether.)

I take my Hermitian inner products to be antilinear in the first entry and linear in the second. If $\mathcal{H}$ is a Hilbert space, the continuous dual of a Hilbert space $\mathcal{H}$ will be denoted $\mathcal{H}^{\vee}$, the complex conjugate space will be denoted $\overline{\mathcal{H}}$. Every linear map $f: \mathcal{H} \rightarrow \mathcal{K}$ is functorially mapped to a linear map $\overline{f}: \overline{\mathcal{H}} \rightarrow \overline{\mathcal{K}}$. If you don't like conjugate spaces, just decorate the correct terms with ``anti".

First note that there is an isometric isomorphism of Hilbert spaces $I: \mathcal{H}_{1} \otimes \mathcal{H}_{2} \rightarrow \mathcal{H}_{1} \otimes \overline{\mathcal{H}_{2}}^{\vee}$ given by linearization and continuous extension of the map $\alpha \otimes \beta \mapsto \alpha \otimes \langle -, \beta \rangle_{2}$. The Hilbert space $\mathcal{H}_{1} \otimes \overline{\mathcal{H}_{2}}^{\vee}$ is canonically identifiable with the space of Hilbert-Schmidt maps from $\overline{\mathcal{H}}_{2}$ to $\mathcal{H}_{1}$; in particular, it is canonically identifiable with a subspace of the space of bounded operators from $\overline{\mathcal{H}}_{2}$ to $\mathcal{H}_{1}$.

Define $X := I(x)$. By the polar decomposition theorem

\begin{align*}
X = U \sqrt{X^{*} X}
\end{align*}

for a partial isometry $U: \overline{\mathcal{H}}_{2} \rightarrow \mathcal{H}_{1}$. The operator $X^{*} X: \overline{\mathcal{H}}_{2} \rightarrow \overline{\mathcal{H}}_{2}$ is positive semidefinite and trace class (with trace $\|x\|^2$), so, it is in particular self-adjoint and compact. By self-adjoint compact operator spectral theory we have that the non-vanishing spectrum $\sigma_{\neq 0}$ of $X^{*} X$ is countable, consists only of eigenvalues, and moreover the eigenspaces $E_{\lambda}$ associated to each $\lambda \in \sigma_{\neq 0}$ are finite dimensional and orthogonal. Positive semidefiniteness ensures that each eigenvalue is positive. By the spectral theorem

\begin{align*}
X^{*} X = \sum_{\lambda \in \sigma_{\neq 0}} \lambda \mathrm{Proj}_{\lambda},
\end{align*}

where $\mathrm{Proj}_{\lambda}$ is the orthogonal projection onto $E_{\lambda}$. Choosing an orthonormal basis $(\beta^{\lambda}_{i})_{i = 1}^{\dim E_{\lambda}}$ for each finite dimensional $E_{\lambda}$

\begin{align*}
X^{*} X &= \sum_{\lambda \in \sigma_{\neq 0}} \lambda \left[\sum_{i_{\lambda} = 1}^{\dim E_{\lambda}} \beta^{\lambda}_{i_{\lambda}} \otimes (\beta^{\lambda}_{i_{\lambda}})^{\vee} \right],
\end{align*}

(where $\omega^{\vee}$ denotes the linear form $\langle \omega, - \rangle$). Of course I can just reindex things here to get a sum

\begin{align*}
X^{*} X=\sum_{j = 1}^{r} \lambda_{j} \beta_{j} \otimes \beta_{j}^{\vee},
\end{align*}

where $r = \mathrm{rank}(X^{*} X) = \sum_{\lambda \in \sigma_{\neq 0}} \dim E_{\lambda} \leq \aleph_{0}$, the $\lambda_{j}$ are eigenvalues, and the $\beta_{j} \in \overline{\mathcal{H}}_{2}$ are associated eigenvectors with $\langle \beta_{i}, \beta_{j} \rangle = \delta_{ij}$. Using our polar decomposition:

\begin{align*}
X = \sum_{j = 1}^{r} \sqrt{\lambda_{j}} (U \beta_{j}) \otimes \beta_{j}^{\vee}.
\end{align*}

Because $U$ is a partial isometry, and one can verify $\ker(U) \leq \ker(X^{*} X)^{\perp}$, the $U \beta_{j} \in \mathcal{H}_{1}$ remain orthonormal. Applying $I^{-1}$, then we simply have
\begin{align*}
x = \sum_{j = 1}^{r} \sqrt{\lambda_{j}} (U \beta_{j}) \otimes \tilde{\beta}_{j}.
\end{align*}
where $\tilde{\beta}_{j} \in \mathcal{H}_{2}$ is the image of $\beta^{\vee}_{j}$ under the canonical isomorphism $\overline{\mathcal{H}}^{\vee}_{2} \overset{\sim}{\rightarrow} \mathcal{H}_{2}$. It is an eigenvector of $\overline{X^{*}X}: \mathcal{H}_{2} \rightarrow \mathcal{H}_{2}$ with eigenvalue $\lambda_{j}$.

A few side remarks:

- $U \beta_{i}$ is an eigenvector of $X X^{*}: \mathcal{H}_{1} \rightarrow \mathcal{H}_{1}$ associated to eigenvalue $\lambda_{i}$. This follows by looking at the explicit construction of $U$ as in, e.g. Douglas' Lemma (a generalization of the polar decomposition theorem). I mention this because $\overline{X^{*} X}$ and $X X^{*}$ are the reduced density states associated to the pure state $x \otimes x^{\vee}$. The proof above is just a generalization of the usual finite-dimensional proofs that the Schmidt decomposition is given by tensor products of eigenstates of reduced density states.
- Because Douglas' lemma constructs a unique $U$, the only choices involved in the procedure above are the choices of bases on each finite dimensional eigenspace $E_{\lambda}$. Any two choices are then related by the action of $\prod_{\lambda \in \sigma_{\neq 0}} U(E_{\lambda})$ (the infinite product of underlying groups). One can check that this group can be realized as a subgroup of $U(\bigoplus_{\lambda \in \sigma_{\neq 0}} E_{\lambda})$, and hence $U(\mathcal{H}_{2})$, or even $U(\mathcal{H}_{1} \otimes \mathcal{H}_{2})$ in a natural way.
- If one wants to be overly pedantic about the question you asked: you mentioned something about bases. Indeed, every Hilbert space has an orthornormal basis, so there exists a (possibly uncountable) basis $(e_{k})_{k \in K}$ for $\mathrm{span}(\beta_{i})^{\perp} \leq \mathcal{H}_{2}$ (I'm taking all spans to be closed) and a basis $(f_{l})_{l \in L}$ for $\mathrm{span}(U\beta_{i})^{\perp} \leq \mathcal{H}_{1}$; we can combine these with the bases we constructed on the spans to get orthonormal bases $(\tilde{e}_{k})_{k \in K'}$ and $(\tilde{f}_{l})_{l \in L'}$. Then, in a rather stupid way, we have
\begin{align*}
x = \sum_{m \in M} \omega_{m} \tilde{e}_{m} \otimes \tilde{f}_{m}
\end{align*}
where $M$ is some (possibly uncountable) indexing set with $|M| = \mathrm{min}(|K'|,|L'|)$. Only countably many of the $\omega_{m}$ are non-zero.