Let $S$ be a stochastic $n$-bit matrix, that is, $S$ is a $(2^n \times 2^n)$-matrix $(a_{ij})$ such that each $a_{ij} \in [0,1]$ and $\sum_{j=0}^{2^n-1} a_{ij} = 1$, for each $i = 0, \dots, 2^n-1$. Writing $|i\rangle$ for the $2^n$-vector representing the binary state of the $n$ bits coding the number $i$, we have $$S|i\rangle = \sum_{j=0}^{2^n-1} a_{ij} |j\rangle,$$ which we interpret as stating that “the probability of outcome $|j\rangle$ upon input $i\rangle$ is $a_{ij}$”.
The goal here is to show how we can simulate the classical probabilistic algorithm coded by $S$ as a quantum algorithm.
The naïve idea is to replace each $a_{ij}$ by $\sqrt{a_{ij}}$; the columns of the resulting matrix $U$ are then all quantum state vectors, and applying the operator $U$ to the basis state $|i\rangle$ produces the quantum state $\sum_{j=0}^{2^n-1} \sqrt{a_{ij}} |j\rangle$. After a von Neumann measurement of this outcome with respect to the computational basis, this means we obtain state $|j\rangle$ with probability $a_{ij}$, which appears to effectively simulate the classical probabilistic algorithm coded by $S$.
However, the problem is that $U$ is not an $n$-qubit operator in general, as it may not be unitary; for instance, the columns of $U$ may not be orthonormal. To simulate the classical probabilistic algorithm coded by $S$ in general, we need “more room” in the form of $n$ ancillary qubits.
So we are looking for a $2n$-qubit operator $U$ that somehow simulates the classical probabilistic algorithm coded by $S$ on the first $n$ qubits.
To do so, we let $|s_i\rangle := \sum_{j=0}^{2^n-1} \sqrt{a_{ij}} |j\rangle$ be the $i$-th column vector of the failed operator $U$ above. Using the Gram-Schmidt procedure, we complete each $|s_i\rangle$ into an orthonormal basis $\{|s_{i0}\rangle = |s_i\rangle, |s_{i1}\rangle, \dots, |s_{i(2^n-1)}\rangle\}$ of the Hilbert space $\CC^{2^n}$. Then we let $U$ be the unique $2n$-qubit operator such that $$U|j\rangle|i\rangle := |i\rangle|s_{ij}\rangle.$$
Remark: In other words, for $i,j \in \{0, \dots, 2^n-1\}$, the $\left(j\cdot2^n + i\right)$-th column of $U$ is $|i\rangle|s_{ij}\rangle$; in particular, the $i$-th column is $|i\rangle|s_i\rangle$.
Let’s check that this $U$ is unitary: given $i,j,k,l \in \{0, \dots, 2^n-1\}$, we have $$\langle i|\langle s_{ij}||k\rangle|s_{kl}\rangle = \langle i|k\rangle \langle s_{ij}|s_{kl}\rangle = \begin{cases} 1 &\text{if } i=k \text{ and } j=l, \\ 0 &\text{otherwise.} \end{cases},$$ as required.
Finally, if we set the ancillary qubits to $|0\rangle$, this $U$ effectively simulates the classical probabilistic algorithm coded by $S$: for $i \in \{0, \dots, 2^n-1\}$, we have $$U|0\rangle|i\rangle = |i\rangle|s_{io}\rangle = |i\rangle|s_i\rangle = \sum_{j=0}^{2^n-1} \sqrt{a_{ij}} |i\rangle|j\rangle.$$ After measuring the first $n$ qubits, this means they are in state $|j\rangle$ with probability $a_{ij}$, as needed.