Re: solving set of inequalities

[Jaroslav Hajek]

On Tue, Jan 26, 2010 at 9:42 AM, Jaroslav Hajek <address@hidden> wrote:
> On Mon, Jan 25, 2010 at 9:41 PM, Petr Korviny <address@hidden> wrote:
>> On 25.1.2010 19:31, John W. Eaton wrote:
>>> On 25-Jan-2010, Petr Korviny wrote:
>>>
>>> | As I mentioned above,
>>>
>>> Above?  Your new text appeared as the first thing in the message I
>>> received, and there was nothing above it.
>>>
>> Sorry, I was writing about text in my earlier posts, my fault.
>>
>>>   I need to find out solution of set of linear inequalities
>>> | (to get values of variables x1,x2,x3,...) or to get an information of
>>> | nonexistence any such solution.
>>> |
>>> | At this picture:
>>> | http://suzelly.opf.slu.cz/~korviny/zzz/octave/excel_solver.png
>>> | you can see utilization of Excel's Solver add-in I used to get solution of
>>> | viewed set. "alpha" variable is always set to a number<0;1>, so set of
>>> | inequalities is linear.
>>>
>>> Then why does your problem statement show you maximizing alpha and
>>> that it has bounds of [0, 1]?
>>>
>>
>> On that picture is statement of complex task and thank you for rewriting 
>> this example
>> with "sqp" function, I certainly use it.
>>
>> In fact I try to divide that problem to two parts/steps:
>> step 1) set variable alpha to constant number (for example to 0.75)
>> step 2) then solve set of linear inequalities, which I get with this alpha
>>
>> According to solution from step 2 I'll change alpha to another number with 
>> new
>> iteration method (step 1) and go to step 2 again. I'll do that until 
>> solution exists for
>> that alpha.
>>
>> My goal is to try to find and optimize this "new iteration method" in step 
>> 1, because the
>> inequalities have some specifics characteristics and we think, me (a 
>> programmer) and my colleague
>> (a mathematician), that the "new iteration method" could be faster than 
>> common used methods.
>>
>> Basic statement is really nonlinear:
>> Max alpha
>> s.t.
>> (1+alpha)*x2 <= x1 <= (3-alpha) * x2
>> (2+alpha)*x3 <= x1 <= (5-2*alpha) * x3
>> 2*x3 <= x2 <= (3-alpha) * x3
>> 0 <= alpha <= 1, x1+x2+x3=1, x1 >= 0, x2 >= 0, x3 >= 0
>>
>> If I set alpha=0 (for simplicity) I get this set of linear inequalities:
>> x2 <= x1
>> x1 <= 3*x2
>> 2*x3 <= x1
>> x1 <= 5*x3
>> 2*x3 <= x2
>> x2 <= 3*x3
>> x1+x2+x3=1, x1 >= 0, x2 >= 0, x3 >= 0
>>
>> And to find some method/function(s), how to solve this with Octave, I'm 
>> trying.
>> I should give you all that information at first probably.
>>
>> Thank you for your time.
>>
>
> I still see one equality. Admittedly, it can be omitted, because all
> inequalities are homogeneous.
> Unless I misunderstood you, however, you don't need to find all
> solutions, just to determine whether a solution exists.
> (I think by default, "solve" means to find all solutions, at least in
> Czech schools :)
> This is obviously a simpler problem:
>
> Let me summarize: you have a matrix
> A = [1, -1, 0; -1, 3, 0; 1, 0, -2; 0, 1, -2; 0, -1, 3];
>
> and you wish to find x so that
> all (A*x >= 0) && all (x >= 0) && sum (x) == 1
> or determine whether it exists.
>
> Preferably using an approach simpler than sqp.
>
> Hmmm.
>
> Perhaps the easiest way is to use random sampling:
>
> [m, n] = size (A);
> x = rand (n, 1000);
> y = A * x;
> solvable = any (all (y >= 0, 1));

Sorry, this should have been:
solvable = any (all (y >= 0, 2));

>
> Something more deterministic that occured to me is to transform it to
> the following quadratic programming problem:
>
> # minimize
> #   norm (y-A*x)^2 + (sum(x) - 1)^2 + eps * norm (x)^2
> # subject to
> #   all (x >= 0) && all (y >= 0)
>
> This is a minimization of a positive definite quadratic form in
> nonnegative variables. For this special kind of quadratic programming
> problems, an efficient algorithm exists that always converges in a
> finite number of steps. It is implemented in Octave 3.3.50+ through
> the function pqpnonneg:
>
> [m, n] = size (A);
> eps = 1e-10; # regularization
> AA = [eye(m), -A; -A', A'*A + ones(n) + eps * eye(n)];
> bb = [zeros(m, 1); -1*ones(n, 1)];
> z = pqpnonneg (AA, bb);
> y = z(1:m);
> x = z(m+1:m+n);
> solvable = norm (A*x - y) <= sqrt(eps);
>
> the disadvantage is that regularization is required to make the
> problem positive definite (it is semidefinite with eps = 0), so even
> if the residual is small, the problem may not be actually solvable,
> only close to being so. Further postprocessing may be required,
> perhaps by sampling in a close neighborhood.
>
> interesting problem. I wish you good luck :)
>
> --
> RNDr. Jaroslav Hajek, PhD
> computing expert & GNU Octave developer
> Aeronautical Research and Test Institute (VZLU)
> Prague, Czech Republic
> url: www.highegg.matfyz.cz
>



-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz