SP heuristic finds incorrect optimum for simple SP

[pgkirsch]

This example should solve to an optimum at x = 1, y = 1.

import gpkit
from gpkit.shortcuts import *

gpkit.enable_signomials() 

x = Var('x')
y = Var('y')

J = 0.01*((x - 1)**2 + (y - 1)**2) + (x*y - 1)**2

sp = SP(J)
print sp
sol = sp.localsolve(x0={x: 1.1, y: 1.1})

print sol.table()

gpkit.disable_signomials()

Output:

In [1]: %run Mike_SP2.py
Optional Python units library (Pint) not installed;unit support disabled.
gpkit.GeometricProgram( # minimize
          -0.02*x + -0.02*y + -2*x*y + 0.01*x**2 + 0.01*y**2 + 1 + x**2*y**2,
          [   # subject to
          ],
          substitutions={  },
          solver="mosek")
Beginning signomial solve.
Using solver 'mosek'
Solving for 2 variables.
Solving took 0.159 seconds.
Solving took 5 GP solves.

Cost
 0.624   

Variables
---------
x : 0.789   
y : 0.789   

Constant sensitivities
----------------------

Possibly related issue: when we print sp it outputs -0.02*x + -0.02*y + -2*x*y + 0.01*x**2 + 0.01*y**2 + 1 + x**2*y**2. Notice that the constant term should be 1.02 instead of 1.

[whoburg]

Todo:

• [ ] fix this
• [ ] add a test
• [ ] also test sol = sp.localsolve() (no initial guess)

[pgkirsch]

Update: The output (specifically the cost) appears to be different from when this ticket was first opened.

[philippekirschen@kirschens-MBP tests]$ ipython issue_274.py 
Optional Python units library (Pint) not installed; unit support disabled.
gpkit.GeometricProgram( # minimize
          -0.02*x + -0.02*y + -2*x*y + 0.01*x**2 + 0.01*y**2 + 1 + x**2*y**2,
          [   # subject to
          ],
          substitutions={  },
          solver="mosek")
Beginning signomial solve.
Using solver 'mosek'
Solving for 2 variables.
Solving took 0.191 seconds.
Solving took 5 GP solves.

Cost
----
 0.1433

Free variables
--------------
x : 0.789 
y : 0.789 

Swept variables
---------------

Constants
---------

Constant and swept variable sensitivities
-----------------------------------------

[pgkirsch]

Observation:

In [28]: gpkit.enable_signomials()

In [29]: x = Var('x')

In [30]: y = Var('y')

In [31]: x-1
Out[31]: gpkit.Signomial(- + x)

In [32]: x-2
Out[32]: gpkit.Signomial(-2 + x)

[pgkirsch]

Another observation:

In [6]: x-0.999
Out[6]: gpkit.Signomial(- + x)

In [7]: x-0.99
Out[7]: gpkit.Signomial(-0.99 + x)

[pgkirsch]

Sub-issues to address:

• [x] The 1 instead of 1.02 constant issue
• [x] The gpkit.Signomial(- + x) issue
• [x] The SP solution issue (if it still remains)

I am open to splitting these into separate issues if anyone has any strong preferences, but for the time being I will keep them here, under the assumption they might be related.

[pgkirsch]

The first sub-issue seems unrelated enough (because it appears to just be printing wrong) and significant enough to merit its own ticket, so I have moved it #316.

[bqpd]

The second issue is happening the same place as the first, I'll move it to #316 as well.

[pgkirsch]

Thanks

[bqpd]

Weirdly, the results on this are better without the 0.01* out front.

[pgkirsch]

And they only get better as the coefficient increases. Fore example, with 10000* in front, the error (cost) is on the order 1e-6.

[pgkirsch]

Getting an error from this

import gpkit
from gpkit.shortcuts import *
gpkit.enable_signomials()

x = Var('x')
y = Var('y')

J = 0.01*((x - 1)**2 + (y - 1)**2) 

sp = SP(J)
sol = sp.localsolve(x0={x: 1.1, y: 1.1})

print sol.table()
gpkit.disable_signomials()

There is a long error traceback sol = sp.localsolve(x0={x: 1.1, y: 1.1}), followed by a rather unsatisfying

KeyError: y

[pgkirsch]

Ok, really weird. I don't get that error without the 0.01 coefficient.

[bqpd]

Not getting that error in the refactor, though it's still getting the wrong answer.

[bqpd]

Hmmm...it gets the correct answer if there's constant term one or greater in the cost. I wonder if the problem is solver limits on driving log(objective) to minus infinity?

[whoburg]

Just a hunch, but I think #296 should take priority -- I wonder if solving that will solve this as well.

[whoburg]

PS confirmed the 1-x issue is fixed on the refactor branch

In [3]: x = Variable('x')

In [4]: from gpkit import enable_signomials

In [5]: enable_signomials()

In [6]: x - 1
Out[6]: gpkit.Signomial(-1 + x)

[pgkirsch]

I will focus on #296 for the time being then.

[pgkirsch]

Unfortunately, commit 35d60b2e635850d6dadadeea984aab7b349b9b33 did not fix this issue too.

[whoburg]

I guess I should have realized it would not -- this fails for both cvxopt and mosek.

[whoburg]

Here's something interesting. Slightly modified version of the original example, with the objective moved to a constraint:

from gpkit import Model, Variable
​
gpkit.enable_signomials() 
x = Variable('x')
y = Variable('y')
z = Variable('z')
​
J = 1*((x - 1)**2 + (y - 1)**2) + (x*y - 1)**2
​
sp = Model(z, [z >= J])
print sp
sol = sp.localsolve()
​
print sol.table()
​
gpkit.disable_signomials()

output:

# minimize
    z,
[   # subject to
    z >= -2*x + -2*x*y + -2*y + 3 + x**2 + x**2*y**2 + y**2,
],
    substitutions={  }
Beginning signomial solve.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-9-e83cd81ff128> in <module>()
     10 sp = Model(z, [z >= J])
     11 print sp
---> 12 sol = sp.localsolve()
     13 
     14 print sol.table()

gpkit/gpkit/model.py in localsolve(self, solver, verbosity, skipfailures, *args, **kwargs)
    279             return self._solve("sp", solver, verbosity, skipfailures, *args, **kwargs)
    280         except ValueError:
--> 281             raise ValueError("'localsolve()' can only be called on models that"
    282                              " contain Signomials, because such"
    283                              " models have only local solutions. Models"

ValueError: 'localsolve()' can only be called on models that contain Signomials, because such models have only local solutions. Models without Signomials have global solutions, so try using 'solve()'.

[bqpd]

Huh, it doesn't for me. It works, actually...

[bqpd]

import gpkit
from gpkit import Model, Variable

gpkit.enable_signomials() 
x = Variable('x')
y = Variable('y')
z = Variable('z')

J = 1*((x - 1)**2 + (y - 1)**2) + (x*y - 1)**2

sp = Model(z, [z >= J])
print sp
sol = sp.localsolve("cvxopt")

print sol.table()

gpkit.disable_signomials()

Outputs:

# minimize
    z,
[   # subject to
    z >= -2*x + -2*x*y + -2*y + 3 + x**2 + x**2*y**2 + y**2,
],
    substitutions={  }
Beginning signomial solve.
Solving took 50 GP solves and 1.15 seconds.

Cost
----
 0.0003349

Free Variables
--------------
x : 0.9999    
y : 0.9999    
z : 0.0003349 

[pgkirsch]

Yah works for me too

[whoburg]

Well that is certainly strange. Here's what's creating the ValueError on my machine (I edited the source code to remove the try/catch in localsolve and raise the error directly):

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-1-75cb0e5eb14d> in <module>()
     11 sp = Model(z, [z >= J])
     12 print sp
---> 13 sol = sp.localsolve("cvxopt")
     14 
     15 print sol.table()

gpkit/gpkit/model.py in localsolve(self, solver, verbosity, skipfailures, *args, **kwargs)
    276         RuntimeWarning if an error occurs in solving or parsing the solution.
    277         """
--> 278         return self._solve("sp", solver, verbosity, skipfailures, *args, **kwargs)
    279             #raise ValueError("'localsolve()' can only be called on models that"
    280             #                  " contain Signomials, because such"

gpkit/gpkit/model.py in _solve(self, programType, solver, verbosity, skipfailures, *args, **kwargs)
    375                 result = self.program.solve(*args, **kwargs)
    376             elif programType == "sp":
--> 377                 result = self.program.localsolve(*args, **kwargs)
    378             sol = self.parse_result(result, constants, mmaps)
    379             solution.append(sol)

gpkit/gpkit/signomial_program.py in localsolve(self, solver, verbosity, x0, reltol, iteration_limit, *args, **kwargs)
    113                and (not (cost and prevcost)
    114                     or abs(prevcost-cost)/(prevcost + cost) > reltol)):
--> 115             gp = self.step(x0, verbosity)
    116             self.gps.append(gp)
    117 

gpkit/gpkit/signomial_program.py in step(self, x0, verbosity)
    165 
    166         gp = GeometricProgram(posy_approxs[0], posy_approxs[1:],
--> 167                               verbosity=verbosity-1)
    168         return gp
    169 

gpkit/gpkit/geometric_program.py in __init__(self, cost, constraints, verbosity)
     44         self.cs = np.hstack((mag(p.cs) for p in self.posynomials))
     45         if not all(self.cs > 0):
---> 46             raise ValueError("GeometricPrograms cannot contain Signomials.")
     47         self.exps = functools_reduce(add, (p.exps for p in self.posynomials))
     48         self.varlocs, self.varkeys = locate_vars(self.exps)

ValueError: GeometricPrograms cannot contain Signomials.

[bqpd]

Hm, looks like a negative c snuck through? Does the code I posted work for you? (you can change the solver, works with all of 'em here)

[whoburg]

Just tried your code above using the latest on master; error is slightly different but still a problem:

python issue274.py
# minimize
    z,
[   # subject to
    z >= -2*x + -2*x*y + -2*y + 3 + x**2 + x**2*y**2 + y**2,
],
    substitutions={  }
Beginning signomial solve.
.../gpkit/gpkit/geometric_program.py:143: RuntimeWarning: overflow encountered in exp
  np.exp(solver_out['primal']).ravel()))
.../gpkit/gpkit/small_scripts.py:38: RuntimeWarning: invalid value encountered in double_scalars
  e = mag(x0[vk]*units * p.diff(vk).sub(x0, require_positive=False).c / p0)
.../gpkit/gpkit/geometric_program.py:51: RuntimeWarning: invalid value encountered in greater
  if not all(self.cs > 0):
Traceback (most recent call last):
  File "issue274.py", line 13, in <module>
    sol = sp.localsolve("cvxopt")
  File ".../gpkit/gpkit/model.py", line 286, in localsolve
    raise ValueError("'localsolve()' can only be called on models that"
ValueError: 'localsolve()' can only be called on models that contain Signomials, because such models have only local solutions. Models without Signomials have global solutions, so try using 'solve()'.

[pgkirsch]

@whoburg are you only getting that error with cvxopt?

[whoburg]

negative, I get it for both cvxopt and mosek. Looks like it's raising early in localsolve before the solver even gets run.

[whoburg]

ps I've removed .pyc files, so it's not that.

[whoburg]

shooting from the hip here, but could be interesting to add this as a unit test on a branch, then we'll get three more machines (on jenkins) trying it out

[bqpd]

Sounds good. I don't think it should be a test long-term, though: takes too long.

[pgkirsch]

Goddamnit @bqpd

[pgkirsch]

This one is at least partly my fault for not saying I was working on it.

[pgkirsch]

Interesting, passed mosek on one machine failed on all other machine/solver combinations.

[bqpd]

Super wacky.

[whoburg]

Agreed re: not making a test; #329 is just debugging

[pgkirsch]

This might be a separate issue, but I thought this was interesting

import gpkit
from gpkit import Model, Variable
with gpkit.EnableSignomials():
    x = Variable('x')
    z = Variable('z')
    J = 0.01*(x - 1)**2
    sp = Model(z, [z >= J])
    sol = sp.localsolve(iteration_limit=50)

#output

Beginning signomial solve.
{x: 0.099503719020998152, z: 0.0094544313448217156}
{x: 0.30858339675657204, z: 0.0070641292715984388}
{x: 0.55138385186126571, z: 0.004137253654887871}
{x: 0.7558555212372341, z: 0.0020329156059758929}
{x: 0.88772160560639424, z: 0.0008742012957197896}
{x: 0.95412334675043098, z: 0.00034540249534560776}
{x: 0.98237655634014498, z: 0.00013081534457576836}
{x: 0.99340766018349569, z: 4.86639535899045e-05}
{x: 0.99755961776297863, z: 1.7977255827085814e-05}
{x: 0.99910015431588806, z: 6.6236742434287643e-06}
{x: 0.99966868273464993, z: 2.4381000785678302e-06}
{x: 0.99987807689494357, z: 8.9711476347932211e-07}
{x: 0.99995514181077461, z: 3.3005551447630321e-07}
{x: 0.9999834968924961, z: 1.2142408125806915e-07}
{x: 0.99999392875103243, z: 4.4669889151738637e-08}
{x: 0.99999776649946515, z: 1.6433196923939157e-08}
{x: 0.99999917833933127, z: 6.0454438375448823e-09}
{x: 0.9999996977276967, z: 2.2239956551815462e-09}
{x: 0.99999988880020207, z: 8.1816243629095626e-10}
{x: 0.99999995909187633, z: 3.0098516157414271e-10}
{x: 0.9999999849507416, z: 1.1072625650178e-10}
{x: 0.99999999446368726, z: 4.0733912983436912e-11}
{x: 1.0, z: 1.0, \fbox{0}: 1.0}
{x: 0.099503719020998152, z: 0.0094544313448217156}
{x: 0.30858339675657204, z: 0.0070641292715984388}
{x: 0.55138385186126571, z: 0.004137253654887871}
{x: 0.7558555212372341, z: 0.0020329156059758929}
{x: 0.88772160560639424, z: 0.0008742012957197896}
{x: 0.95412334675043098, z: 0.00034540249534560776}
{x: 0.98237655634014498, z: 0.00013081534457576836}
{x: 0.99340766018349569, z: 4.86639535899045e-05}
{x: 0.99755961776297863, z: 1.7977255827085814e-05}
{x: 0.99910015431588806, z: 6.6236742434287643e-06}
{x: 0.99966868273464993, z: 2.4381000785678302e-06}
{x: 0.99987807689494357, z: 8.9711476347932211e-07}
{x: 0.99995514181077461, z: 3.3005551447630321e-07}
{x: 0.9999834968924961, z: 1.2142408125806915e-07}
{x: 0.99999392875103243, z: 4.4669889151738637e-08}
{x: 0.99999776649946515, z: 1.6433196923939157e-08}
{x: 0.99999917833933127, z: 6.0454438375448823e-09}
{x: 0.9999996977276967, z: 2.2239956551815462e-09}
{x: 0.99999988880020207, z: 8.1816243629095626e-10}
{x: 0.99999995909187633, z: 3.0098516157414271e-10}
{x: 0.9999999849507416, z: 1.1072625650178e-10}
{x: 0.99999999446368726, z: 4.0733912983436912e-11}
{x: 1.0, z: 1.0, \fbox{1}: 1.0}
{x: 0.099503719020998152, z: 0.0094544313448217156}
{x: 0.30858339675657204, z: 0.0070641292715984388}
{x: 0.55138385186126571, z: 0.004137253654887871}
{x: 0.7558555212372341, z: 0.0020329156059758929}
Solving took 50 GP solves and 0.828 seconds.

Cost
----
 0.002033

Free Variables
--------------
x : 0.7559   
z : 0.002033 

Unsurprisingly, if you set iteration_limit=22, the answer is correct:

Cost
----
 4.073e-11

Free Variables
--------------
x : 1         
z : 4.073e-11 

tl;dr the while loop exit condition in localsolve doesn't seem to be doing its job.

[bqpd]

There's a cycle! Haha, weird...maybe some of those are feasibility solves? (which aren't permitted to be the last solve)

[pgkirsch]

Alas no it doesn't seem like the except statement in localsolve ever gets activated so they are all plain vanilla gp.solves.

[bqpd]

Weirdddd. Is the cost difference really staying above reltol?

[pgkirsch]

Is this what you mean? (the try is necessary for when prevcost=None)

 try:
       print abs(prevcost-cost)/(prevcost + cost)
except:
       pass

in which case yes...

# output

0.144704016817
0.261295916399
0.341050295372
0.398578505591
0.433582450504
0.450607123003
0.457720705538
0.460476303345
0.461510257991
0.461893445602
0.462034797781
0.462086850697
0.462106007011
0.462113054959
0.462115648176
0.462116602485
0.462116949706
0.462117078157
0.462117132357
0.462117195347

[pgkirsch]

By the way please can you explain the first part of that conditional?

        while (iterations < iteration_limit
               and (not (cost and prevcost)
                    or abs(prevcost-cost)/(prevcost + cost) > reltol)):

(starting from not)

[bqpd]

Oh, that's catching when those costs are None, such as the first two iterations and whenever a feasibility gp is used.

Re: costs, hmmmmmm. Perhaps we should make that denominator (prevcost + cost + 1e-10?)

[pgkirsch]

Doesn't seem to help unfortunately

[bqpd]

A large enough denominator certainly would :p I think the problem is that prevcost+cost is approaching 0 just as fast as cost is approaching 0

[pgkirsch]

Oh wait you were right, I think the end of each "cycle" is a feasible solve...

[pgkirsch]

Putting 1e-6 in the denominator does indeed solve the problem.

Moving this discussion to a separate issue...

[pgkirsch]

@whoburg @bqpd any suggestions for things I can try, I'm running out of ideas...

[bqpd]

The feasible-solve is interesting; what's the gp just before that in model.program? (it should be the infeasible one)

[pgkirsch]

Mm no sorry, I meant the original issue, I think you're referring to what is now #336.