I need to generate a list of random variables $\bf{x}$ subject to constraints that can be expressed in the form $\bf{E}x=b$ where $\bf{E}$ is an $m \times n $ matrix if $\bf{x}$ has $n$ entries. In all the cases I'm dealing with, $n >> m$, for example $n$ will be around 14,000 and $m$ will be 50. I'm not sure what method I will use for random sampling, either normal or uniform, it is not clear which is best for the problem I'm trying to solve, but I need all the variables to be sampled from distributions with the same mean and range/variance.
What I've been doing to solve this is reducing $\bf{E}$ to row-echelon form, setting all the variables corresponding to columns to the right of the last pivot to random values, and then solving the remaining square matrix equality.
There is a problem however, to solve the square matrix equality, I subtract the already set values from the right hand side. Unfortunately, the variances add as well, so my last 50 values tend to vary hugely, which is unfortunately unacceptable in this problem.
Is there a better way to do this? I cannot think of a way to fix the current method I am using. I use R.