Consider an arbitrary matrix $M \in \mathbb{R}^{n \times m}$. Denote the column space of $M$ as $\mathcal{C}(M)$.
Is it always possible to construct a right stochastic matrix $S$ such that $\mathcal{C}(M) \subseteq \mathcal{C}(S)$?
Of course the new matrix $S$ can have more columns. Ideally, I would like a procedure for constructing $S$ that uses a reasonably small number of columns (i.e. linear in $m$).
Yes. If $m > n/2$, then $n$ is $O(m)$, and you can take the matrix $S$ to be the identity.
If $m < n/2$, then you can do the following:
The columns of $M$ form a spanning set for $C(M)$, right? Let's call the $i$th column $c_i$. Split $c_i = p_i + n_i$, where $p_i$ has all the positive-or-zero entries, and $n_i$ has the negative entries. (I.e., $p_i$ is $c_i$ with all negative entries replaced by zeroes). Divide each $p_i$ and $n_i$ by the absolute value of its largest entry, and then divide each by $2m$; negate the $n_i$ vector as well. Call the resulting vectors $p'_i$ and $n'_i$ Now the matrix $H$, whose columns are all the $p'_i$ and $n'_i$ vectors, has the same column space as $M$, and $2n < m$ columns. But it's not right-stochastic: many rows may sum to a number less than one. (Each entry in each column is less than $m/2$). Let $v$ be the sum of the columns, and let $S$ be the matrix $H$ with the column $1 - v$ appended, and then enough zero-columns to make it square. Then $C(S)$ contains $C(M)$, and has $2m+1 = O(m)$ nonzero columns.
Example: $$ M = \begin{bmatrix} 3 & -1 \\ 1 & 1 \\ -4 & 0 \\ 2 & 0 \\ 0 & 1 \end{bmatrix} $$ Then $m = 2, n = 3$. And initially, \begin{align} p_1 = \begin{bmatrix} 3 \\ 1 \\ 0 \\ 2 \\ 0 \end{bmatrix}, n_1 = \begin{bmatrix} 0 \\ 0 \\ -4 \\ 0\\ 0 \end{bmatrix}, p_2 = \begin{bmatrix} 0 \\ 1 \\ 0 \\ 0 \\ 1 \end{bmatrix}, n_2 = \begin{bmatrix} -1 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \end{align} After division by the largest absolute entry of each vector, and by $2m = 4$, and negation of the $n_i$s, we have \begin{align} p'_1 = \begin{bmatrix} 1/4 \\ 1/12 \\ 0 \\ 1/6\\ 0 \end{bmatrix}, n'_1 = \begin{bmatrix} 0 \\ 0 \\ 1/4 \\ 0\\ 0 \end{bmatrix}, p'_2 = \begin{bmatrix} 0 \\ 1/4 \\ 0\\ 0\\ 1/4 \end{bmatrix}, n'_2 = \begin{bmatrix} 1/4 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \end{align} The matrix $H$ is then \begin{align} \begin{bmatrix} 1/4 & 0 & 0 & 1/4 \\ 1/12 & 0 & 1/4 & 0 \\ 0 & 1/4 & 0 & 0 \\ 1/6 & 0 & 0 & 0 \\ 0 & 0 & 1/4 & 0 \end{bmatrix} \end{align} The column sum for $H$ is the vector \begin{align} \begin{bmatrix} 1/2 \\ 1/3 \\ 1/4 \\ 1/6 \\ 1/4 \end{bmatrix} \end{align} and when we append this to $H$ we get the matrix $S$ (with no zero-columns, because $n = 2m+1$, rather than being larger).
$$ S = \begin{bmatrix} 1/4 & 0 & 0 & 1/4 & 1 - 1/2\\ 1/12 & 0 & 1/4 & 0 & 1 - 1/3\\ 0 & 1/4 & 0 & 0 & 1 - 1/4\\ 1/6 & 0 & 0 & 0 & 1 - 1/6\\ 0 & 0 & 1/4 & 0 & 1 - 1/4 \end{bmatrix} $$
It's true that in this case, $S$ might as well be the identity, but I wanted to show you the process for a real example, and this was the first that came to mind.