Let $X$ be $n$-dimensional, $Y$ be $m$-dimensional, and both with zero means. (It's simple to deal with nonzero means later because they just get added in.) We seek an $n$-dimensional zero-mean random vector $Z$ with unit covariance matrix and uncorrelated with $Y$ (that is, $E[ZY^\prime]=0$), an $n\times m$ matrix $A,$ and an $n\times n$ matrix $B$ for which
$$X = AY + BZ.$$
Rewriting this as
$$X - AY = BZ,$$
right-multiplying by the row vector $Y^\prime$ and taking expectations gives
$$\operatorname{Cov}(X,Y) - A\operatorname{Cov}(Y) = E[XY^\prime] - AE[YY^\prime] = BE[ZY^\prime] = 0.\tag{*}$$
When $\operatorname{Cov}(Y)$ is nonsingular this has a unique solution
$$A = \operatorname{Cov}(X,Y)\operatorname{Cov}(Y)^{-1}.$$
Take covariances of the $n$-vectors in the rewritten equation to find
$$\operatorname{Cov}(X) - A\operatorname{Cov}(X,Y) - \operatorname{Cov}(Y,X)A^\prime + A\operatorname{Cov}(Y)A^\prime = BB^\prime.$$
This exhibits $B$ as a(ny) matrix square root of the (symmetric, positive semidefinite) matrix on the left hand side.
This solution assumed nothing about the actual distribution of $(X,Y,Z);$ but when it is multivariate Normal, the vanishing covariance of $Z$ with $Y$ assures $Z$ and $Y$ are also independent and because $AY+BZ$ is $n$-variate Normal and has the same second moments as $X,$ it gives the intended distribution.
When the covariance of $Y$ is singular, there might still be a solution. In such circumstances you will have to examine $(*)$ more closely, for instance by row-reducing the $m\times(n+m)$ augmented matrix $[\operatorname{Cov}(Y) \mid \operatorname{Cov}(Y,X)].$