0
$\begingroup$

Link : Reduce time taken to get feasible solution

I have 18 inequations with 4 variables vars2 = {arb[1], arb[2], arb[3], arb[4]} , To get a feasible solution, NMinimize and Nmaximize were used as was suggested in the link. But unfortunately this method is not returning any feasible solution (which was required).

expr2 = (31 (52275 arb[1] + 17425 arb[2] - 
   8659 (3 arb[3] + arb[4])))/(
181178952 (-arb[2] arb[3] + arb[1] arb[4])) >= 0 && (
31 (6531 arb[1] + 2177 arb[2] - 4283 (3 arb[3] + arb[4])))/(
181178952 (-arb[2] arb[3] + arb[1] arb[4])) >= 0 && (
31 (129 arb[1] + 43 arb[2] - 13 (3 arb[3] + arb[4])))/(
4213464 (-arb[2] arb[3] + arb[1] arb[4])) >= 0 && (
505325 arb[1] + 5035825 arb[2] - 8659 (29 arb[3] + 289 arb[4]))/(
60392984 (arb[2] arb[3] - arb[1] arb[4])) >= 0 && (
63133 arb[1] + 629153 arb[2] - 4283 (29 arb[3] + 289 arb[4]))/(
60392984 (arb[2] arb[3] - arb[1] arb[4])) >= 
0 && (-1247 arb[1] - 12427 arb[2] + 377 arb[3] + 3757 arb[4])/(
1404488 (-arb[2] arb[3] + arb[1] arb[4])) >= 
0 && (-557600 arb[1] + 2247825 arb[2] + 
 8659 (32 arb[3] - 129 arb[4]))/(
90589476 (arb[2] arb[3] - arb[1] arb[4])) >= 
0 && (-69664 arb[1] + 280833 arb[2] + 
 4283 (32 arb[3] - 129 arb[4]))/(
90589476 (arb[2] arb[3] - arb[1] arb[4])) >= 0 && (
1376 arb[1] - 5547 arb[2] - 416 arb[3] + 1677 arb[4])/(
2106732 (-arb[2] arb[3] + arb[1] arb[4])) >= 
0 && ((93 - 745 arb[2]) arb[3] + (31 + 745 arb[1]) arb[4])/(
838 (-arb[2] arb[3] + arb[1] arb[4])) >= 0 && (
31 (-arb[2] (-1 + arb[3]) + arb[1] (3 + arb[4])))/(
838 (arb[2] arb[3] - arb[1] arb[4])) >= 0 && (
31 (arb[2] (29 - 774 arb[3]) + 31 (3 arb[3] + arb[4]) + 
   arb[1] (87 + 774 arb[4])))/(
216204 (-arb[2] arb[3] + arb[1] arb[4])) >= 0 && (
87 (1 + arb[2]) arb[3] + 3 (289 - 29 arb[1]) arb[4])/(
838 (arb[2] arb[3] - arb[1] arb[4])) >= 
0 && (-867 arb[2] (-1 + arb[3]) + arb[1] (87 + 867 arb[4]))/(
838 (-arb[2] arb[3] + arb[1] arb[4])) >= 0 && (
arb[2] (8381 - 7482 arb[3]) + 899 arb[3] + 8959 arb[4] + 
 29 arb[1] (29 + 258 arb[4]))/(
72068 (arb[2] arb[3] - arb[1] arb[4])) >= 
0 && ((32 - 387 arb[2]) arb[3] + 129 (-1 + 3 arb[1]) arb[4])/(
419 (-arb[2] arb[3] + arb[1] arb[4])) >= 0 && (
129 arb[2] (-1 + arb[3]) + arb[1] (32 - 129 arb[4]))/(
419 (arb[2] arb[3] - arb[1] arb[4])) >= 0 && (
992 arb[3] - 129 arb[2] (29 + 64 arb[3]) - 3999 arb[4] + 
 32 arb[1] (29 + 258 arb[4]))/(
108102 (-arb[2] arb[3] + arb[1] arb[4])) >= 0;

vars2 = arb /@ Range[4];

cons2 = And @@ Thread[-10 <= vars2 <= 10];

(soln3 = NMinimize[{Total[List @@ expr2[[All, 1]]], expr2 && cons2}, 
 vars2, Reals]) // Quiet // AbsoluteTiming

expr2 && cons2 /. soln3[[2]]

(soln4 = NMaximize[{Total[List @@ expr2[[All, 1]]], expr2 && cons2}, 
 vars2, Reals]) // Quiet // AbsoluteTiming

expr2 && cons2 /. soln4[[2]]

To make sure there exist some feasible solution, I have examples which were found using randomSearch technique (took a long time).

{arb[1] -> 151/2000, arb[2] -> 88147/258000, arb[3] -> 4019/4000, arb[4] -> -(67417/516000)}

{arb[1] -> 49/160, arb[2] -> 46837/64500, arb[3] -> 7071/4000, arb[4] -> 25619/129000}

soleg1=Thread[vars2 -> {49/160, 46837/64500, 7071/4000, 25619/129000}];
soleg2=Thread[vars2 -> {151/2000, 88147/258000, 4019/4000, -(67417/516000)}];
expr2 && cons2 /. soleg1 (* Returns True *)
expr2 && cons2 /. soleg2 (* Returns True *)

As asked in the comments about groupInverse of a matrix A. Code to find group Inverse of a matrix of index 1 :

function Ahash = phash( A )
%Matrix return group inverse of A
%Date:1/10/2017

  r=rank(A);
  if(r==rank(A*A))
    Q=[orth(A) null(A,'r')];
    P=Q\A*Q;
    C=P(1:r,1:r);
    D=zeros(size(A,1));
    D(1:r,1:r)=inv(C);
    inv(C)
    Ahash=Q*D/Q;
  else
    disp('Invalid input : Matrix is not of index 1.');
    Ahash=-1;
  end
end
Improve this question
$\endgroup$
9
  • $\begingroup$ So, what is the question? $\endgroup$ – Henrik Schumacher Oct 19 '17 at 9:19
  • $\begingroup$ ... But unfortunately this method (as suggested in link) is not returning any feasible solution $\endgroup$ – Pushpendu Ghosh Oct 19 '17 at 9:56
  • $\begingroup$ The question is I want to find feasible solutions of expr2, in an optimized faster way (RandomSearch takes lot of time) $\endgroup$ – Pushpendu Ghosh Oct 19 '17 at 9:58
  • $\begingroup$ Where do these inequalities come from? Seemingly from 2 x 2 matrices, since the denominators (-arb[2] arb[3] + arb[1] arb[4]) look like determinants. I am not sure if this helps: If you add the requirement that you only want to study matrices with nonnegative determinant, i.e., by add (-arb[2] arb[3] + arb[1] arb[4]) >=0 to the constraints and use Simplify, you get much simpler polynomial inequalities. And {0,0,0,0} is a feasible point for that. $\endgroup$ – Henrik Schumacher Oct 19 '17 at 10:51
  • $\begingroup$ What is my main aim? Given any singular matrix A(nxn) , st groupInverse(A) is non negative. I need to find a matrix U (called Proper splitting of A) st. 1. Range(A) = Range(U) and 2. NullSpace(A)=NullSpace(U). Here 1. and 2. give linear relations between the elements of U, solving it I get some independent variables which here, in this code is vars2. Hence if I know numerical values of vars2, I get a U which is a proper splitting of A. Now out of these conditions - top half are : groupInverse(U) is non negative (n^2 elements) and rest are : groupInverse(U).(U-A) is non negative (n^2 elements). $\endgroup$ – Pushpendu Ghosh Oct 19 '17 at 11:24
1
$\begingroup$

Since your equations are greaterequal zero, for a solution they must have a definite values b1,b2,...

The trick is to bring the inequations into a linear form where there are no quotients. This is much easier to handle for Mathematica. Therefore work with equations as follows: 

list5 = {(
31 (52275 arb[1] + 17425 arb[2] - 8659 (3 arb[3] + arb[4])))/(
181178952 (-arb[2] arb[3] + arb[1] arb[4])) == b1, (
31 (6531 arb[1] + 2177 arb[2] - 4283 (3 arb[3] + arb[4])))/(
181178952 (-arb[2] arb[3] + arb[1] arb[4])) == b2, (
31 (129 arb[1] + 43 arb[2] - 13 (3 arb[3] + arb[4])))/(
4213464 (-arb[2] arb[3] + arb[1] arb[4])) == b3, (
505325 arb[1] + 5035825 arb[2] - 8659 (29 arb[3] + 289 arb[4]))/(
60392984 (arb[2] arb[3] - arb[1] arb[4])) == b4, (
63133 arb[1] + 629153 arb[2] - 4283 (29 arb[3] + 289 arb[4]))/(
60392984 (arb[2] arb[3] - arb[1] arb[4])) == 
b5, (-1247 arb[1] - 12427 arb[2] + 377 arb[3] + 3757 arb[4])/(
1404488 (-arb[2] arb[3] + arb[1] arb[4])) == 
b6, (-557600 arb[1] + 2247825 arb[2] + 
 8659 (32 arb[3] - 129 arb[4]))/(
90589476 (arb[2] arb[3] - arb[1] arb[4])) == 
b7, (-69664 arb[1] + 280833 arb[2] + 
 4283 (32 arb[3] - 129 arb[4]))/(
90589476 (arb[2] arb[3] - arb[1] arb[4])) == b8, (
1376 arb[1] - 5547 arb[2] - 416 arb[3] + 1677 arb[4])/(
2106732 (-arb[2] arb[3] + arb[1] arb[4])) == 
b9, ((93 - 745 arb[2]) arb[3] + (31 + 745 arb[1]) arb[4])/(
838 (-arb[2] arb[3] + arb[1] arb[4])) == b10, (
31 (-arb[2] (-1 + arb[3]) + arb[1] (3 + arb[4])))/(
838 (arb[2] arb[3] - arb[1] arb[4])) == b11, (
31 (arb[2] (29 - 774 arb[3]) + 31 (3 arb[3] + arb[4]) + 
   arb[1] (87 + 774 arb[4])))/(
216204 (-arb[2] arb[3] + arb[1] arb[4])) == b12, (
87 (1 + arb[2]) arb[3] + 3 (289 - 29 arb[1]) arb[4])/(
838 (arb[2] arb[3] - arb[1] arb[4])) == 
b13, (-867 arb[2] (-1 + arb[3]) + arb[1] (87 + 867 arb[4]))/(
838 (-arb[2] arb[3] + arb[1] arb[4])) == b14, (
arb[2] (8381 - 7482 arb[3]) + 899 arb[3] + 8959 arb[4] + 
 29 arb[1] (29 + 258 arb[4]))/(
72068 (arb[2] arb[3] - arb[1] arb[4])) == 
b15, ((32 - 387 arb[2]) arb[3] + 129 (-1 + 3 arb[1]) arb[4])/(
419 (-arb[2] arb[3] + arb[1] arb[4])) == b16, (
129 arb[2] (-1 + arb[3]) + arb[1] (32 - 129 arb[4]))/(
419 (arb[2] arb[3] - arb[1] arb[4])) == b17, (
992 arb[3] - 129 arb[2] (29 + 64 arb[3]) - 3999 arb[4] + 
 32 arb[1] (29 + 258 arb[4]))/(
108102 (-arb[2] arb[3] + arb[1] arb[4])) == b18};

Eliminate the arb[i]

eli = Eliminate[list5, {arb[1], arb[2], arb[3], arb[4]}]

(*     93 b1 == -29 b5 - 258 b6 + 87 b8 + 774 b9 && 
  77934 b10 == 
  69285 + 1549462 b5 - 1824326 b6 - 4648386 b8 + 5472978 b9 && 
  77934 b11 == -2883 - 468442 b5 + 3589154 b6 + 1405326 b8 - 
  10767462 b9 && 
  77934 b12 == 8649 + 238830 b5 - 622634 b6 - 716490 b8 + 1867902 b9 &&
  25978 b13 == 2697 - 1549462 b5 + 1824326 b6 && 
  25978 b14 == 26877 + 468442 b5 - 3589154 b6 && 
  25978 b15 == -2697 - 238830 b5 + 622634 b6 && 
  12989 b16 == 11997 - 774731 b8 + 912163 b9 && 
  12989 b17 == 3999 + 234221 b8 - 1794577 b9 && 
  12989 b18 == 992 - 119415 b8 + 311317 b9 && 3 b2 == -b5 + 3 b8 && 
  3 b3 == -b6 + 3 b9 && 31 b4 == 29 b5 + 258 b6 && 
  31 b7 == 29 b8 + 258 b9     *)

Solve for the bi. You have four free parameters, b3, b7, b11, b17.

{sol = Solve[eli]} // TableForm

(*     {{b13 -> 109875510/47165573 + (129 b11)/13 - (52275 b17)/8659 +    (
15084 b3)/13 - (648612 b7)/8659, 
  b5 -> -(136591725/4056239278) - (93 b11)/559 + (774 b17)/8659 - (
12849 b3)/559 + (12849 b7)/8659, 
  b6 -> 11223/3628121 - (87 b17)/8659 - 3 b3 + (1677 b7)/8659, 
  b15 -> 567690135/2028119639 + (855 b11)/559 - (9201 b17)/8659 + (
77934 b3)/559 - (77934 b7)/8659, 
  b4 -> -(2697/468442) - (87 b11)/559 - (25977 b3)/559 + 3 b7, 
  b14 -> -3 b11 + 3 b17, 
  b2 -> 961/468442 + (31 b11)/559 + (4283 b3)/559, 
  b8 -> -(33282/3628121) + (258 b17)/8659 + (4283 b7)/8659, 
  b10 -> 4176/5447 - (43 b11)/13 - (5028 b3)/13, 
  b9 -> 3741/3628121 - (29 b17)/8659 + (559 b7)/8659, 
  b12 -> 21576/234221 - (285 b11)/559 - (25978 b3)/559, 
  b16 -> 5598858/3628121 - (17425 b17)/8659 - (216204 b7)/8659, 
  b18 -> 672731/3628121 - (3067 b17)/8659 - (25978 b7)/8659, 
  b1 -> 899/468442 + (29 b11)/559 + (8659 b3)/559}}     *)

With this bi solve for the arb[i]

sol2 = Solve[list5 /. First@sol, {arb[1], arb[2], arb[3], arb[4]}]

(*      {{arb[1] -> -(112567 (258 b11 + 31 b17))/(12 (-115971 - 3134958 b11 + 
     376681 b17 + 10182538 b11 b17 - 936055218 b3 + 
     3040365398 b17 b3 - 7260851 b7 - 196277198 b11 b7)), 
  arb[2] -> -(112567 (31 + 64 b11 - 93 b17))/(12 (-115971 - 
     3134958 b11 + 376681 b17 + 10182538 b11 b17 - 936055218 b3 + 
     3040365398 b17 b3 - 7260851 b7 - 196277198 b11 b7)), 
  arb[3] -> (-1391652 - 96062946 b11 - 7022275 b17 - 11232662616 b3 - 
  87130212 b7)/(12 (-115971 - 3134958 b11 + 376681 b17 + 
    10182538 b11 b17 - 936055218 b3 + 3040365398 b17 b3 - 
    7260851 b7 - 196277198 b11 b7)), 
  arb[4] -> (-7367491 - 23829568 b11 + 21066825 b17 - 2786396928 b3 + 
  261390636 b7)/(12 (-115971 - 3134958 b11 + 376681 b17 + 
    10182538 b11 b17 - 936055218 b3 + 3040365398 b17 b3 - 
    7260851 b7 - 196277198 b11 b7))}}     *)

Inserting these arb[i] into the original equations helps to FindInstance of the bi quite fast and the corresponding arb[i]

fi11 = FindInstance[
   expr2 /. First@sol2 // Simplify, {b3, b7, b11, b17}, 2]

(*     {{b3 -> 57916370079/283236901694603, b7 -> 787/23911, 
   b11 -> 53401995667357991381/392803698272339835314, 
   b17 -> 269/1419}, {b3 -> 2110071648/16201787509361, b7 -> 61/2403, 
   b11 -> 6527487720/38667750619, b17 -> 952600891/3088030419}}     *)

sol11 = First@sol2 /. First[fi11]

(*     {arb[1] -> 1181512951903198958775209827554373/
    4303233111631058231748194542658916, 
       arb[2] -> 106128231536901263581238652240929/
    717205518605176371958032423776486, 
       arb[3] -> 5368850238953974904842591393832689/
    4303233111631058231748194542658916, 
      arb[4] -> -(30330695488083688169589434961005/
    358602759302588185979016211888243)}     *)

Test the original equations

expr2 /. sol11 // Simplify

(+     True     *)

Now, with Reduce you get conditions for the bi much faster (not shown here) as if reducing the original equations for the arb[i].

Reduce[expr2 /. First@sol2 // Simplify, {b3, b7, b11, b17}]

Get an impression, how the arb[i] behave:

(fi12 = FindInstance[
 expr2 /. First@sol2 // Simplify, {b3, b7, b11, b17}, 
 100]); // Timing       (* about 100 seconds   *)

sol12 = First@sol2 /. fi12;

ListLinePlot[{arb[1], arb[2], arb[3], arb[4]} /. sol12]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.