Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Speeding-up FindInstance

• To: mathgroup at smc.vnet.net
• Subject: [mg94954] Re: [mg94919] Speeding-up FindInstance
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 31 Dec 2008 06:11:02 -0500 (EST)
• References: <200812301055.FAA13175@smc.vnet.net>

Colletti, Bruce W. wrote:
> Re Mathematica 7.0.
> 
> Let X be an n-by-p array of indexed x-variables, e.g., x[1,1], x[1,2], etc.
> 
> I seek x-values in {-1, 0, 1} so that all column sums are zero and no column is the zero vector.
> 
> The code below for a 5x3 problem runs forever.  What is a faster way to get a solution?  Thankx.
> 
> Bruce
> 
> {n,p} = {5,3};
> 
> X=Array[x,{n,p}];
> 
> Z = Flatten[X];
> FindInstance[And@@Join[Map[#==0&,Total@X],
> Map[#>0&,Map[Abs,Total@X,{2}]],
> Map[Abs@#<=1&,Z]],Z,Integers]

I trust this is a simplified example of a harder class of problem? For 
this specific setup you could just take a row of -1's, a row of 1's, and 
all other rows just 0's.

Anyway, you'll want to make it explicitly a linear problem. This means 
the absolute values need to go away. You can do this be defining a 
second set of "shadow" variables, and imposing constraints that force 
them to be the absolute values of the desired set. The code below will 
do this, with a slight catch. It remains not quite an integer linear 
program due to use of inequations (linear forms set not_equal to one). 
There might be a better way to do it but I've not stumbled across it as 
yet (still working on coffee #2).

{n, p} = {5, 3};
X = Array[x, {n, p}];
absX = Array[ax, {n, p}];
Z = Flatten[X];
aZ = Flatten[absX];

In[40]:= constraints =
  Join[Map[# == 0 &, Total@X], Map[# >= 1 &, Total@absX],
   Thread[aZ >= Z], Thread[aZ >= -Z], Map[-1 <= # <= 1 &, Z],
   Thread[aZ + Z != 1], Map[0 <= # <= 1 &, aZ]]

Out[40]= {x[1, 1] + x[2, 1] + x[3, 1] + x[4, 1] + x[5, 1] == 0,
  x[1, 2] + x[2, 2] + x[3, 2] + x[4, 2] + x[5, 2] == 0,
  x[1, 3] + x[2, 3] + x[3, 3] + x[4, 3] + x[5, 3] == 0,
  ax[1, 1] + ax[2, 1] + ax[3, 1] + ax[4, 1] + ax[5, 1] >= 1,
  ax[1, 2] + ax[2, 2] + ax[3, 2] + ax[4, 2] + ax[5, 2] >= 1,
  ax[1, 3] + ax[2, 3] + ax[3, 3] + ax[4, 3] + ax[5, 3] >= 1,
  ax[1, 1] >= x[1, 1], ax[1, 2] >= x[1, 2], ax[1, 3] >= x[1, 3],
  ax[2, 1] >= x[2, 1], ax[2, 2] >= x[2, 2], ax[2, 3] >= x[2, 3],
  ax[3, 1] >= x[3, 1], ax[3, 2] >= x[3, 2], ax[3, 3] >= x[3, 3],
  ax[4, 1] >= x[4, 1], ax[4, 2] >= x[4, 2], ax[4, 3] >= x[4, 3],
  ax[5, 1] >= x[5, 1], ax[5, 2] >= x[5, 2], ax[5, 3] >= x[5, 3],
  ax[1, 1] >= -x[1, 1], ax[1, 2] >= -x[1, 2], ax[1, 3] >= -x[1, 3],
  ax[2, 1] >= -x[2, 1], ax[2, 2] >= -x[2, 2], ax[2, 3] >= -x[2, 3],
  ax[3, 1] >= -x[3, 1], ax[3, 2] >= -x[3, 2], ax[3, 3] >= -x[3, 3],
  ax[4, 1] >= -x[4, 1], ax[4, 2] >= -x[4, 2], ax[4, 3] >= -x[4, 3],
  ax[5, 1] >= -x[5, 1], ax[5, 2] >= -x[5, 2],
  ax[5, 3] >= -x[5, 3], -1 <= x[1, 1] <= 1, -1 <= x[1, 2] <= 1, -1 <=
   x[1, 3] <= 1, -1 <= x[2, 1] <= 1, -1 <= x[2, 2] <= 1, -1 <=
   x[2, 3] <= 1, -1 <= x[3, 1] <= 1, -1 <= x[3, 2] <= 1, -1 <=
   x[3, 3] <= 1, -1 <= x[4, 1] <= 1, -1 <= x[4, 2] <= 1, -1 <=
   x[4, 3] <= 1, -1 <= x[5, 1] <= 1, -1 <= x[5, 2] <= 1, -1 <=
   x[5, 3] <= 1, ax[1, 1] + x[1, 1] != 1, ax[1, 2] + x[1, 2] != 1,
  ax[1, 3] + x[1, 3] != 1, ax[2, 1] + x[2, 1] != 1,
  ax[2, 2] + x[2, 2] != 1, ax[2, 3] + x[2, 3] != 1,
  ax[3, 1] + x[3, 1] != 1, ax[3, 2] + x[3, 2] != 1,
  ax[3, 3] + x[3, 3] != 1, ax[4, 1] + x[4, 1] != 1,
  ax[4, 2] + x[4, 2] != 1, ax[4, 3] + x[4, 3] != 1,
  ax[5, 1] + x[5, 1] != 1, ax[5, 2] + x[5, 2] != 1,
  ax[5, 3] + x[5, 3] != 1, 0 <= ax[1, 1] <= 1, 0 <= ax[1, 2] <= 1,
  0 <= ax[1, 3] <= 1, 0 <= ax[2, 1] <= 1, 0 <= ax[2, 2] <= 1,
  0 <= ax[2, 3] <= 1, 0 <= ax[3, 1] <= 1, 0 <= ax[3, 2] <= 1,
  0 <= ax[3, 3] <= 1, 0 <= ax[4, 1] <= 1, 0 <= ax[4, 2] <= 1,
  0 <= ax[4, 3] <= 1, 0 <= ax[5, 1] <= 1, 0 <= ax[5, 2] <= 1,
  0 <= ax[5, 3] <= 1}

In[41]:= Timing[
  result = FindInstance[constraints, Join[Z, aZ], Integers]]

Out[41]= {24.8062, {{x[1, 1] -> 1, x[1, 2] -> 1, x[1, 3] -> 1,
    x[2, 1] -> 1, x[2, 2] -> 1, x[2, 3] -> 1, x[3, 1] -> 0,
    x[3, 2] -> 0, x[3, 3] -> -1, x[4, 1] -> -1, x[4, 2] -> -1,
    x[4, 3] -> -1, x[5, 1] -> -1, x[5, 2] -> -1, x[5, 3] -> 0,
    ax[1, 1] -> 1, ax[1, 2] -> 1, ax[1, 3] -> 1, ax[2, 1] -> 1,
    ax[2, 2] -> 1, ax[2, 3] -> 1, ax[3, 1] -> 0, ax[3, 2] -> 0,
    ax[3, 3] -> 1, ax[4, 1] -> 1, ax[4, 2] -> 1, ax[4, 3] -> 1,
    ax[5, 1] -> 1, ax[5, 2] -> 1, ax[5, 3] -> 0}}}

In[42]:= X /. result
Out[42]= {{{1, 1, 1}, {1, 1, 1}, {0, 0, -1}, {-1, -1, -1}, {-1, -1,
    0}}}

An alternative is to make it a quadratic program, using squares instead 
of emulating absolute values. For this particular example the timing 
seems about the same. I have no idea what to expect in general.

In[47]:= constraints2 =
  Join[Map[# == 0 &, Total@X], Thread[Total[X^2] >= 1],
   Map[-1 <= # <= 1 &, Z]]

Out[47]= {x[1, 1] + x[2, 1] + x[3, 1] + x[4, 1] + x[5, 1] == 0,
  x[1, 2] + x[2, 2] + x[3, 2] + x[4, 2] + x[5, 2] == 0,
  x[1, 3] + x[2, 3] + x[3, 3] + x[4, 3] + x[5, 3] == 0,
  x[1, 1]^2 + x[2, 1]^2 + x[3, 1]^2 + x[4, 1]^2 + x[5, 1]^2 >= 1,
  x[1, 2]^2 + x[2, 2]^2 + x[3, 2]^2 + x[4, 2]^2 + x[5, 2]^2 >= 1,
  x[1, 3]^2 + x[2, 3]^2 + x[3, 3]^2 + x[4, 3]^2 + x[5, 3]^2 >= 1, -1 <=
    x[1, 1] <= 1, -1 <= x[1, 2] <= 1, -1 <= x[1, 3] <= 1, -1 <=
   x[2, 1] <= 1, -1 <= x[2, 2] <= 1, -1 <= x[2, 3] <= 1, -1 <=
   x[3, 1] <= 1, -1 <= x[3, 2] <= 1, -1 <= x[3, 3] <= 1, -1 <=
   x[4, 1] <= 1, -1 <= x[4, 2] <= 1, -1 <= x[4, 3] <= 1, -1 <=
   x[5, 1] <= 1, -1 <= x[5, 2] <= 1, -1 <= x[5, 3] <= 1}

In[49]:= Timing[result2 = FindInstance[constraints2, Z, Integers]]

Out[49]= {23.4654, {{x[1, 1] -> -1, x[1, 2] -> -1, x[1, 3] -> -1,
    x[2, 1] -> -1, x[2, 2] -> -1, x[2, 3] -> -1, x[3, 1] -> 0,
    x[3, 2] -> 0, x[3, 3] -> 0, x[4, 1] -> 1, x[4, 2] -> 1,
    x[4, 3] -> 1, x[5, 1] -> 1, x[5, 2] -> 1, x[5, 3] -> 1}}}

Daniel Lichtblau
Wolfram Research

• References:
  • Speeding-up FindInstance
    • From: "Colletti, Bruce W." <bcolletti@nvcc.edu>