Overview
The short answer is no you cannot prepare this state efficiently if you wish to know the values of $y_x$; but, in fact, your question is a special case of a more general case I prove below which I will outline first as it is applicable to a more general audience. After I will show that your question is indeed a special case of the theorem. Finally, I will address the comments under your original post concerning mixed states and how in some cases the use of a mixed state may allow for efficient preparation with the caveat that the values of $y_x$ will not be known.
Theorem
Let $\mathcal H_\alpha$ and $\mathcal H_\beta$ be the Hilbert spaces containing states with subscripts $\alpha$ and $\beta$ respectively where $N\equiv\operatorname{dim}\mathcal H_\alpha$. Suppose we wish to prepare the state:
$$|\psi\rangle=\sum_{x\in X}a_x|x\rangle_\alpha|z_x\rangle_\beta$$
where $\mathcal B\equiv\left\{|x\rangle\middle|\forall x\in X\right\}$ is a complete orthonormal basis of $\mathcal H_\alpha$; and $a_x\ne0$ for all $x$. Let $\mathcal Z\equiv\left(|z_x\rangle\middle|\forall x\in X\right)$ be an $N$-tuple that contains all duplicates not just unique states in an order defined by $x$. Suppose further, we already have the multiset of states $\mathcal S\equiv\left\{|0\rangle_\alpha\right\}\cup\left\{|z_x\rangle\in\mathcal Z\right\}$ in quantum memory i.e. we have 1 copy of $\mathcal H_\alpha$ and $N$ copies of $\mathcal H_\beta$. Let $\mathcal A$ be some set of $\mathcal Z$ such that knowledge of a single element from $\mathcal Z\in\mathcal A$, it's position in $\mathcal Z$ and knowledge of $\mathcal A$ gives no information on rest of the states in $\mathcal Z$. i.e. $\mathcal A$ can be written as the cartesian product of sets $\mathcal D_i\subseteq\mathcal H_\beta$ each with cardinality greater than one:
$$\mathcal A\equiv\prod_{i=1}^N\mathcal D_i\quad\text{where }\left|\mathcal D_i\right|>1$$
It is impossible to prepare $|\psi\rangle$ from the states in $\mathcal S$ in better than $\mathcal O\left(N\right)$ complexity in 1,2-qubit gates for a given $\mathcal A$ using the same quantum circuit for all $\mathcal Z\in\mathcal A$.
Proof
To prepare $|\psi\rangle$ one must transfer the information from the states in $\mathcal Z$ into our single state in $\mathcal H_\alpha$. Each $|z_x\rangle$ is independent - given one $|z_x\rangle$ we have no information about the rest of $\mathcal Z$. Due to this independence clearly any operation involving only one $|z_x\rangle$ can only transfer the information contained within that specific $|z_x\rangle$ and as any quantum operation can be broken down into one and 2-qubit gates; then we can consider the transfer operation for each $|z_x\rangle$ as a separate operation each with one or more 2-qubit gates. Note one or more 2-qubit gates because a 1-qubit gate clearly cannot transfer information from one qubit to another; however, a 2-qubit gate can transfer the information from one qubit to another (consider the SWAP gate). But as there are $N$ states in $\mathcal Z$ then we require at least $N$ 2-qubit gates.
Preperation of State in Question
Now, quantum circuits can comprise three types of actions: addition of ancilla lines, unitary operations and measurements. Ancilla lines can always be added at the beginning of the circuit and in and of themselves will not introduce any randomness allowing for the production of $\left\{y_x\right\}$. unitary operations are deterministic and so cannot introduce randomness either. Thus, the randomness must come from measurements.
To use measurements to sample from a probability distribution we must prepare a state that when we measure it we obtain each $y_x$ with its corresponding probability from $D$ - or at least measure some state that can be mapped onto $y_x$ uniquely. If you want $N$ independent values of $y_x$ then you must sample $D$ at least $N$ times. Alternatively, we could sample the distribution $D$ without the use of the quantum computer $N$ times and then prepare the $N$ states with this knowledge.
Now we sample $N$ times to get the $N$-tuple of states $\left(|y_1\rangle,\ldots,|y_N\rangle\right)$ or at least $\left(\left|\gamma_1\right\rangle,\ldots,\left|\gamma_N\right\rangle\right)$ where a unitary $U$ maps from $|\gamma\rangle\to|y\rangle$. If it was more efficient to prepare the state that gives $\left\{|\gamma\rangle\right\}$ with probability distribution $D$ than $\left\{|y\rangle\right\}$ then it will be more efficient to apply $U$ once after we prepare $\frac{1}{\sqrt{N}}\sum_x|x\rangle|\gamma_x\rangle$ rather than $N$ times before. But either way we now want to prepare a state of the form given in the above theorem where $\mathcal D_i=\left\{|\gamma\rangle\right\}$ where in the case $U=I$ then $|\gamma\rangle=|y\rangle$. Thus, by the above theorem it must take at least $N$ 2-qubit gates to prepare $\frac{1}{\sqrt{N}}\sum_x|x\rangle|\gamma_x\rangle$. Therefore, there is no efficient way to produce $|\psi\rangle$ if $D$ has more than one non-zero probability. Additionally, to this it may also be inefficient to produce the initial state from which you measure to produce the $N$-tuple $\left(\left|\gamma_1\right\rangle,\ldots,\left|\gamma_N\right\rangle\right)$.
Note I have neither specified the dimensions of the Hilbert space within which $|y\rangle$ reside; whether $|y\rangle$ are product states, entangled states or even within the computational basis; nor whether $D$ is efficiently sampleable or not in order to make the proof more general.
Pure or Mixed
A mixed state encompasses our classical uncertainty of which state the system is in. A mixed state $\rho$ is given by:
$$\rho\equiv\sum_iP_i|\psi_i\rangle\langle\psi_i|$$
where $P_i$ is the probability of $|\psi_i\rangle$ being the quantum state of the system. Thus the action of a unitary operator on a mixed state is given by $\rho'=U^\dagger\rho U$. Now, let $U|\psi_i\rangle=|\phi_i\rangle$ for all $i$. If we substitute this in we find that:
$$\rho'=\sum_iP_i|\phi_i\rangle\langle\phi_i|$$
So the action of the unitary $U$ on a mixed state is the same as the action of $U$ on each pure state comprising $\rho$ and then combining these pure states with the same probabilities.
Thus, whether we use a mixed or pure state to describe the ensemble of possible $|\psi\rangle$ is really immaterial to the action of the circuit and so using the proof above one still cannot prepare the desired state efficiently if one wishes to know the values of $y_x$ drawn from $D$. However, the decomposition of mixed states as a sum of pure states is not unique and this allows us to increase the efficiency in some cases.
If you wish to represent $|\psi\rangle$ as a mixed state then you would write:
$$\rho=\sum_Yp\left(Y\right)|\psi_Y\rangle\langle\psi_Y|$$
where $Y\equiv\left(y_1,\ldots,y_N\right)$ and $p\left(Y\right)\equiv\prod_{i=1}^Np\left(y_i\right)$ is the probability of getting $Y$ and $p\left(y_i\right)$ is the probability in $D$ of getting $y_i$.
The mixed state $\rho$ has $E\ge N\left(|D|-1\right)$ non-zero eigenvalues (lower bound proved below) and so it is possible to express $\rho$ as:
$$\rho=\sum_{i=1}^E\rho_i|\rho_i\rangle\langle\rho_i|$$
where $\left\{|\rho_i\rangle\right\}$ are mutually orthonormal eigenvectors of $\rho$ with corresponding non-zero eigenvalue $\rho_i$ which in principle can be calculated.
Now if we can implement a unitary $U_\rho$ that that rotates the computational basis into the eigenbasis of $\rho$ (and its inverse) efficiently, and additionally can prepare the state $|\phi\rangle:=\sum_{i=1}^E\sqrt{\rho_i}|\rho_i\rangle$ efficiently then the following process will prepare the mixed state $\rho$ efficiently:
- prepare $|\phi\rangle$
- apply $U_\rho^{-1}$
- measure in the computational basis
- finally, apply $U_\rho$ to the collapsed state to efficiently produce $\rho$
This is the same as measuring $|\phi\rangle$ in the eigenbasis of $\rho$. However, this will in general not be efficient and in this case, we do not know then set $\left\{y_x\right\}$ of values drawn from $D$.
Proof of lower bound on the number of eigenvalues
The number of eigenvalues may be of help in determining the complexity and happened to be something I proved while thinking about the problem so I have included it below:
Consider:
$$\left(\langle z|\langle r|\right)|y\rangle=\frac{1}{\sqrt{N}}\sum_{x\in X}\langle z|x\rangle\langle r|y_x\rangle=\frac{1}{\sqrt{N}}\langle r|y_z\rangle$$
Thus, the matrix elements of $\rho$ are:
$$\begin{align}\langle z|\langle r|\rho|x\rangle|s\rangle&=\frac{1}{N}\sum_{y_1\in D}\ldots\sum_{y_N\in D}p\left(y_1\right)\ldots p\left(y_N\right)\langle r|y_z\rangle\langle y_x|s\rangle\\&=\frac{1}{N}\begin{cases}p_rp_s,&x\ne z\\p_r,&x=z\text{ and }r=s\\0,&x=z\text{ and }r\ne s\end{cases}\end{align}$$
where $p_y:=p\left(y\right)$ for ease of notation.
Therefore, we can write $\rho$ as a block matrix with the diagonal matrix $A:=\frac{1}{N}\operatorname{diag}\left(p_1,\ldots,p_{|D|}\right)$ along the diagonal of $\rho$ and the projector $B:=\frac{1}{N}|p\rangle\langle p|$ on the off diagonals where $|p\rangle:=\sum_{i\in D}p_i|i\rangle$ and $|D|$ is the number of values that $y_x$ can take on with non-zero probability ($|D|=N$ in the question):
$$\rho=\begin{bmatrix}A&B&B&\dots&B\\B&A&B&\dots&B\\B&B&A&\dots&B\\\vdots&\vdots&\vdots&\ddots&\vdots\\B&B&B&\dots&A\end{bmatrix}$$
Consider the null space of $\rho$:
$$\begin{align}\rho\begin{bmatrix}|x_1\rangle\\\vdots\\|x_N\rangle\end{bmatrix}&=0\\\iff A|x_i\rangle+B\sum_{\substack{j=1\\i\ne j}}^N|x_j\rangle&=0\quad\forall i\in\left[1,N\right]\\\iff A|x_i\rangle&=-\frac{1}{N}|p\rangle\langle p|\sum_{\substack{j=1\\i\ne j}}^N|x_j\rangle\quad\forall i\in\left[1,N\right]\end{align}$$
As $A^{-1}=N\operatorname{diag}\left(\frac{1}{p_1},\ldots,\frac{1}{p_{|D|}}\right)$ then $|x_i\rangle\propto\sum_{i\in D}|i\rangle$. This leaves $N$ degrees of freedom, the constants of proportionality, so the dimensionality of the null space of $\rho$ is at most $N$ and so $E\ge N\left(|D|-1\right)$.