# Multiple solutions for linear equations

Can I generate multiple (random) vectors [x1, x2, x3, x4, x5] where all x's are positive numbers such that they all satisfy the linear constraints:

x1+2.x2+6.x3+4.x4+2.x5=10

3.x1+x2+4.x3+7.x4+4.x5=16

Can I for example generate 100 instances of [x1, x2, x3, x4, x5]? And how?

Call $$\mathfrak A$$ the set $$\mathfrak A=\left\{x\in\mathbb R_+^5;Ax=b \right\}$$ where $$A=\left(\begin{matrix} 1 &2 &6 &4 &2\\ 3 &1 &4 &7 &4\end{matrix} \right)\qquad b={10 \choose 16}$$ It is the intersection of two hyperplanes of $$\mathbb R^5$$ restricted to the positive quadrant. One can then generate a Markov chain on $$\mathfrak A$$ as follows:

1. Start from an arbitrary $$x^{(0)}\in\mathfrak A$$
2. For iteration $$t$$, elect at random three different integers in $$\{1,2,3,4,5\}$$, $$i_1,i_2,i_3$$
3. Generate $$y_{i_1},y_{i_2},y_{i_3}$$ uniformly in $$[0,10]^3$$
4. If there exist a solution to $$Ay=b$$ with the three above components, switch $$x^{(t)}$$ to this solution $$y^*$$, else repeat $$x^{(t)}$$
5. Increment $$t$$ and go to 2.

This Markov chain should have uniform distribution over $$\mathfrak A$$.

Here is for instance an R code implementation:

#warning: codegolfed!
m=rbind(m<-1:5,o<-m)
A=rbind(c(1,2,6,4,2),c(3,1,4,7,4))
b=c(10,16)
u=c(6,5,2,3,4)
for(t in 1:1e4){
o[i<-sample(1:5,3)]=y=runif(3,0,u[i])
p=ifelse(rep(!prod((o[-i]<-solve(B<-A[,-i],
b-B%*%y))),5),m[length(m[,1]),],o)
m=rbind(m,p)}


An alternative is to consider the uniform distribution over $$\begin{cases} x_1+2x_2+6x_3\le 10\\ 3x_1+x_2+4x_3\le16\end{cases}\tag{1}$$ which can be simulated by drawing points uniformly in $$[0,16/3]\times[0,5]\times[0,5/3]$$ and keeping only those satisfying (1).

There are two equations with five unknowns so that in general you can solve only two unknowns. For example, you can solve $$x_4$$ and $$x_5$$ from the system. Now you only need to sample $$x_1$$, $$x_2$$, and $$x_3$$. Since they are positive, you can sample from the range as given by @Xi'an: $$[0, 16/3] \times [0,5] \times [0, 5/3]$$. After that you calculate $$x_4$$ and $$x_5$$ and check whether they are positive. If yes, keep that vector, otherwise, sample and calculate again.

EDIT: I have edit the half-close interval into close interval. As required, we have to check whether the sampled values are 0 or positive.

• Yes, what I mean is to check whether they are positive or not. If 0, then the procedure has to be repeated. Dec 10, 2020 at 21:33
• You are correct at least in R as R says "runif will not generate either of the extreme values unless max = min or max-min is small compared to min, and in particular not for the default arguments". But is it correct in general, and how do we know that? I mean, mathematically, (0, 1) is not the same as [0,1]. Do I miss something from Uniform(a, b) or Uniform[a, b]. Could you help to clarify? Dec 11, 2020 at 9:38
• As I stated in the first comment, a Uniform$(a,b)$ is the same thing as a Uniform$([a,b]$ since changing the density over a set of measure zero does not modify the associated distribution. Dec 11, 2020 at 10:44
• Yes, you are right mathematically. But what will be wrong if for example R sets to sample 0 or 1 from U(0,1)? I do not see any disapproval for R to sample the end points? In other words, I do not see why R has to restrict in (0, 1). Of course, it is not mathematically necessary, but a machine can do, right? Dec 11, 2020 at 22:04
• This discussion is quite tangential to the purpose of the original question: please check this thread on stackoverflow. Dec 12, 2020 at 9:02