Finding A Random Solution With All Values Within A Distribution That Fit A Set Of Linear Equations With Less Equations Than Unknowns I have a number of linear equations with more variables than equations, and I am trying to find random solutions to these equations so that all of the values lie within a range. I already know that these values exist and just want to find them. A simpler explanation would be preferable along with pseudocode. Thank you!
(edit) I am asking this question based on neural networks and I am trying to map the outputs of the network to its inputs, and the layer I am inverting has more inputs than outputs, so I only need something that is inputted that outputs the output I am trying to invert. I tried adding equations to the set by setting the weights and biases to be a random distribution of the current weights and biases.
 A: Just to make sure, what you are trying to do: 
You are given a system of equations given by 
$$ Ax=b$$ with $A\in \mathbb{R}^{n\times m}$ and $b\in \mathbb{R}^{n}$ where $m>n$, i.e. an underdetermined system of equations with a wide matrix (more columns than rows). 
You are interested in a "random" solution $x$ with "values in a range". I take "random" to mean arbitrary, i.e. any solution is valid and the range to be elementwise. 
First of all: rewrite your system of equation by shifting and scaling to make the range of values equal for all entries of $x$ and centered around zero. Shifting by a transformation of $b$, scaling by a transformation in $A$ (I can give more details, if required). 
You are then given a system of equations with 
$Ax=b$ with $\|x\|_\infty <q $, where $A$ and $b$ are your transformed matrix and right-hand side and $q$ is your allowed range. 
I propose using something I call the "iterative Euclidian minimization". You start by finding a solution to your linear system with $\|x\|_2$ minimal. This can be done with the pseudo-inverse of $A$ obtained via (for example) SVD. 
This gives a first approximation but some values might lie out of range (i.e $\|x\|_\infty>q$). 
Let $x_i=q_1$ with $|q_1|>q$ be the most extreme of these outliers (i.e the entry with the largest absolute value in your current approximation). You can decrease the norm of this vector by forcing this entry a little bit towards zero and treating this entry as fixed. You add the respective column of $A$ over to $b$ and solve the new linear system with one degree of freedom less (again with the Pseudo-inverse of the new matrix with one column less). This gives a solution to $Ax=b$ with a smaller $\|x\|_\infty$. 
Repeating this procedure of choosing the maximum element (or multiple maximum elements)and shifting them slightly towards zero will decrease the maximum norm of $x$ and will eventually lead to a vector with $\|x\|_\infty<q$, if one exists. 
There are some pitfalls and the question of "how much do I shift them inwards", but I leave that to you to figure it out. 
Update after a comment from OP: 
1) My range is two-sided as well. I set bounds for maximum upper and lower bound for any entry of $x$ in the same manner. Just keep it symmetric around $0$ and watch out for the signs of the entries. 
2) I could have made this clearer here: Once you fixed an entry and obtained another solution, you shift that index further inward (as it is still one of the outliers)
3) This was one of the pitfalls, I was talking about: Shift this new outlier also into this range and fix it. Let the rest of the unfixed variable do the rest
