# Get a 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

• So, what is the question? – Henrik Schumacher Oct 19 '17 at 9:19
• ... But unfortunately this method (as suggested in link) is not returning any feasible solution – Pushpendu Ghosh Oct 19 '17 at 9:56
• The question is I want to find feasible solutions of expr2, in an optimized faster way (RandomSearch takes lot of time) – Pushpendu Ghosh Oct 19 '17 at 9:58
• 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. – Henrik Schumacher Oct 19 '17 at 10:51
• 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). – Pushpendu Ghosh Oct 19 '17 at 11:24

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]