best way to obtain one answer that satisfy a linear equation in matlab
Asked Answered
Q

1

3

I have a linear equation:

vt = v1*x1 + v2*x2 + v3*x3

vt, v1, v2, v3 are scalars with values between 0 and 1. What is the best way to generate one set (any set will be fine) of x1, x2 and x3 that satisfy the equation above. and also satisfy

x1>0
x2>0
x3>0

I have couple thousand sets of vt,v1,v2 and v3, therefore I need to be able to generate x1, x2 and x3 programmatically.

Quinquefid answered 16/7, 2014 at 15:40 Comment(4)
The fastest way would be to set x1 and x2 to be known and solve for x3. x1 and x2 could be set by you, or they could be randomly generated. If you randomly generate them, you should make sure that vt is less than v1*x1 + v2*x2. You could also use linprog and form this into a linear program and solve for it. Does it matter if x1 and x2 are the same?Knott
@Knott This was my initial thought as well. I can't think of a somewhat clever way to randomly generate x1 and x2 and ensure vt<v1*x1 + v2*x2 (I know I can keep randomly generate x1 and x2 until the condition is met).Quinquefid
@Knott also, it is preferable if x1 and x2 are drawn from some distribution (e.g. a uniform distribution) randomly (randomly with some constraint).Quinquefid
That's ok. We can avoid this all together by using linprog. I'll write an answer for you.Knott
K
3

There are two ways you could approach this:

  1. From the method that you have devised in your post. Randomly generate x1 and x2 and ensure that vt < v1*x1 + v2*x2, then go ahead and solve for x3.
  2. Formulate this into linear program. A linear program is basically solving a system of equations that are subject to inequality or equality constraints. In other words:

blah

As such, we can translate your problem to be of a linear programming problem. The "maximize" statement is what is known as the objective function - the overall goal of what you are trying to accomplish. In linear programming problems, we are trying to minimize or maximize this objective. To do this, we must satisfy the inequalities seen in the subject to condition. Usually, this program is represented in canonical form, and so the constraints on each variable should be positive.

The maximize condition can be arbitrary as you don't care about the objective. You just care about any solution. This whole paradigm can be achieved by linprog in MATLAB. What you should be careful with is how linprog is specified. In fact, the objective is minimized instead of maximized. The conditions, however, are the same with the exception of ensuring that all of the variables are positive. We will have to code that in ourselves.

In terms of the arbitrary objective, we can simply do x1 + x2 + x3. As such, c = [1 1 1]. Our equality constraint is: v1*x1 + v2*x2 + v3*x3 = vt. We also must make sure that x is positive. In order to code this in, what we can do is choose a small constant so that all values of x are greater than this value. Right now, linprog does not support strict inequalities (i.e. x > 0) and so we have to circumvent this by doing this trick. Also, to ensure that the values are positive, linprog assumes that the Ax <= b. Therefore, a common trick that is used is to negate the inequality of x >= 0, and so this is equivalent to -x <= 0. To ensure the values are non-zero, we would actually do: -x <= -eps, where eps is a small constant. However, when I was doing experiments, by doing it this way, two of the variables end up to be the same solution. As such, what I would recommend we do is to generate good solutions that are random each time, let's draw b to be from a uniform random distribution as you said. This will then give us a starting point every time we want to solve this problem.

Therefore, our inequalities are:

 -x1 <= -rand1
 -x2 <= -rand2
 -x3 <= -rand3

rand1, rand2, rand3 are three randomly generated numbers that are between 0 and 1. In matrix form, this is:

 [-1 0 0][x1]      [-rand1]
 [0 -1 0][x2]  <=  [-rand2]
 [0 0 -1][x3]      [-rand3]

Finally, our equality constraint from before is:

[v1 v2 v3][x1]     [vt] 
          [x2]  = 
          [x3]

Now, to use linprog, you would do this:

X = linprog(c, A, b, Aeq, beq);

c is a coefficient array that is defined for the objective. In this case, it would be defined as [1 1 1], A and b is the matrix and column vector defined for the inequality constraints and Aeq and beq is the matrix and column vector defined for the equality constraints. X would thus give us the solution after linprog converges (i.e. x1, x2, x3). As such, you would do this:

A = -eye(3,3);
b = -rand(3,1);
Aeq = [v1 v2 v3];
beq = vt;
c = [1 1 1];
X = linprog(c, A, b, Aeq, beq);

As an example, suppose v1 = 0.33, v2 = 0.5, v3 = 0.2, and vt = 2.5. Therefore:

rng(123); %// Set seed for reproducibility
v1 = 0.33; v2 = 0.5; v3 = 0.2;
vt = 2.5;
A = -eye(3,3);
b = -rand(3,1);
Aeq = [v1 v2 v3];
beq = vt;
c = [1 1 1];
X = linprog(c, A, b, Aeq, beq);

I get:

X =

0.6964
4.4495
0.2268

To verify that this equals vt, we would do:

s = Aeq*X

s = 2.5000

The above simply does v1*x1 + v2*x2 + v3*x3. This is computed in a dot product form to make things easy as X is a column vector and v1, v2, v3 are already set in Aeq and is a row vector.


As such, either way is good, but at least with linprog, you don't have to keep looping until you get that condition to be satisfied!

Small Caveat

One small caveat that I forgot to mention in the above approach is that you need to make sure that vt >= v1*rand1 + v2*rand2 + v3*rand3 to ensure convergence. Since you said that v1,v2,v3 are bounded between 0 and 1, the worst case is when v1,v2,v3 are all equal to 1. As such, we really need to make sure that vt > rand1 + rand2 + rand3. If this is not the case, then simply take each value of rand1, rand2, rand3, and divide by (rand1 + rand2 + rand3) / vt. As such, this will ensure that the total summation will equal vt assuming that all of the weights are 1, and this will allow the linear program to converge properly.

If you don't, then the solution will not converge due to the inequality conditions placed in for b, and you won't get the right answer. Just some food for thought! As such, do this for b before you run linprog

if sum(-b) > vt
   b = b ./ (sum(-b) / vt);
end

Good luck!

Knott answered 16/7, 2014 at 16:35 Comment(6)
Thanks for the answer! One question, so does vector b corresponds to the -eps? in that case does it make sense to do something like b = -rand(3,1)*0.0001? @KnottQuinquefid
@Quinquefid - Yes, which is what I noted in my Small Caveat. Multiply the rand by something small to ensure that we get an answer. Other than that, this should hopefully give you what you're looking for.Knott
now I see why b = -rand(3,1)*0.0001 is not a super great idea! since sometime it just return an answer that is too close to zero.Quinquefid
@Quinquefid - Which is why you should scale it according to what I have specified in my Small Caveat. That way, you won't get solutions so close to 0.Knott
Thanks for the great answer! I should add that, your answer is super easy to extend to a system that have N number of variables in the linear equation. super awesome! @KnottQuinquefid
@Quinquefid - You're very welcome! I made a small correction to the Small Caveat. I didn't divide properly. Good luck!Knott

© 2022 - 2024 — McMap. All rights reserved.