Topic: Unexpected inequality constraints in Ipopt with FuncDesigner

I'm experimenting with solving Thomson problem (electrons with Coulomb potential on a sphere) from http://www.mcs.anl.gov/~more/cops/ using OpenOpt with Ipopt solver. However I get very bad performance when I try to use AD. Here is the code:

from openopt import NLP
from numpy import *
from FuncDesigner import *
from math import pi

np = 50 # number of electrons

# Set the starting point to a quasi-uniform distribution
# of electrons on a unit sphere

random.seed(861276191)
theta = 2*pi*random.random_sample(np)
phi = pi*random.random_sample(np)

x, y, z = oovar(size=np), oovar(size=np), oovar(size=np)

f = sum([sum([1.0/sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2 + (z[i] - z[j])**2) \
    for j in range(i+1, np)]) for i in range(np-1)])

x0 = {x : cos(theta)*sin(phi), y : sin(theta)*sin(phi), z : cos(phi)}

p = NLP(f, x0)
p.constraints = [x**2 + y**2 + z**2 == ones(np)]

r = p.solve('ipopt', options = 'output_file = result.txt, file_print_level = 5')
x_opt, y_opt, z_opt = r(x,y,z)
print(x_opt, y_opt, z_opt)

Ipopt output:

...

Number of nonzeros in equality constraint Jacobian...:     7500
Number of nonzeros in inequality constraint Jacobian.:   183750
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:      150
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       50
Total number of inequality constraints...............:     1225
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:     1225

...

Total CPU secs in IPOPT (w/o function evaluations)   =    985.358
Total CPU secs in NLP function evaluations           =   1562.662

I can't understand why 1225 inequality constraints appear. I suppose they slow down the solution process greatly. Another strange thing is 7500 nonzeros in Jacobian of equality constraints. Because there should be only 50 nonzeros. Are this two anomalies ok? Or maybe I've used OpenOpt incorrectly?

I use OpenOpt 0.32 with Python 2.6.6 and Ipopt 3.8.3 with sequential MUMPS 4.9.2 on Debian Squeeze.

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

http://openopt.org/FuncDesignerDoc#Attached_constraints

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

Thank you. I've tried adding removeAttachedConstraints() call but Ipopt still outputs the same number of constraints and spends as much time as it was.

...

f = sum([sum([1.0/sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2 + (z[i] - z[j])**2) \
    for j in range(i+1, np)]) for i in range(np-1)])

f.removeAttachedConstraints()

x0 = {x : cos(theta)*sin(phi), y : sin(theta)*sin(phi), z : cos(phi)}

p = NLP(f, x0)
p.constraints = [x**2 + y**2 + z**2 == ones(np)]

...

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

you should remove from each sqrt personally, e.g.
S = []
for i in xrange(N):
    s = []
    for j in xrange(M):
        tmp_x,tmp_y,tmp_z = sqrt(...), sqrt(...),sqrt(...)
        tmp_x.removeAttachedConstraints()
        tmp_y.removeAttachedConstraints()
        tmp_z.removeAttachedConstraints()
        s.append(1/tmp_x + 1/tmp_y + 1/tmp_z)
    S += sum(s)
f = sum(S)
Pay attention that your unvectorized problem will work slow in Python interpreter cycles. Maybe involving Unladen Swallow could speedup it.

Regards, D.

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

Another approach (quite ugly although):
p=NLP(...)
p._prepare()
p.c = None
r=p.minimize(...)

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

here's my a little bit vectorized version:

from openopt import NLP
from numpy import *
from FuncDesigner import *
from math import pi

np = 50 # number of electrons

# Set the starting point to a quasi-uniform distribution
# of electrons on a unit sphere

random.seed(861276191)
theta = 2*pi*random.random_sample(np)
phi = pi*random.random_sample(np)

points = oovars(np) # np points of size 3 for each one
S = []
for i in range(np-1):
    for j in range(i+1, np) :
        tmp = sqrt(sum((points[i]-points[j])**2))
        tmp.removeAttachedConstraints()
        S.append(1/tmp)
f = sum(S)

points_start = dict([(points[i], (cos(theta[i])*sin(phi[i]), sin(theta[i])*sin(phi[i]), cos(phi[i]))) for i in xrange(np)])
p = NLP(f, points_start)
p.constraints = [sum(points[i]**2) == 1 for i in xrange(np)]
r = p.solve('ipopt', useSparse=False, options = 'output_file = result.txt, file_print_level = 5')

I used ipopt 3.8.3; I guess 3.9.0 would work better, but I have some troubles in linking it with pyipopt.
So, my results (AMD Athlon 3800):
iter    objFunVal    log10(maxResidual)   
    0  1.493e+03             -15.65
   10  1.071e+03              -0.08
   20  1.056e+03              -2.73
   30  1.055e+03              -3.44
   40  1.055e+03              -3.35
   50  1.055e+03              -4.68
   60  1.055e+03              -4.97
   70  1.055e+03              -6.00
   80  1.055e+03              -9.60
   90  1.055e+03              -9.13
   97  1.055e+03             -10.90
istop: 1000
Solver:   Time Elapsed = 292.45     CPU Time Elapsed = 270.59
objFunValue: 1055.1823 (feasible, MaxResidual = 1.26452e-11)
however, using random numbers yields different elapsed time.
BTW, you could import pi from numpy, not Python math module.

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

Very impressive! After doing removeAttachedConstraints correctly it started to work 5.5 times faster. And 2.3 times faster after vectorization. 12x speedup. Great. Thank you!

I didn't know that it's possible to make a vectorized solution. It looks much more natural (exclusive of dict usage). Actually I've barely debugged the original one (due to == operator in constraints which was hard to find). However it was half a year ago.

Another interesting thing is that in the vectorized program we've got

Number of nonzeros in equality constraint Jacobian...:      150

That is zeros were found correctly by the OpenOpt. Which was not the case for the previous solution. Why?

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

As for nonzeros in constraint Jacobian, unfortunately OpenOpt currently can't construct dependency matrix correctly for vectorized function.
E.g. for
a=oovar(size=50)
c=a**2<1
openopt supplies dependency matrix (called by IPOPT before optimization process start) will be full of ones instead of being diagonal. It can be replaced by
a=oovars(50)
and sometimes this is helpful for low-RAM computers, but sometimes it will work slower, and sometimes the vectorized operations like the ones I used here will not work.
Generally it is recommended to use oovar() instead of oovars(), and RAM issue is the only one recommendation to use oovars() instead.

Here's even more vectorized approach, it creates only O(n) auxiliary oofuns to be summarized instead of O(n^2)
By the way, I wrote incorrectly: it is xrange to be removed in Python3, and range will return iterator instead of a list.

from openopt import NLP
from numpy import *
from FuncDesigner import *
from math import pi

np = 50 # number of electrons

# Set the starting point to a quasi-uniform distribution
# of electrons on a unit sphere

random.seed(861276191)
theta = 2*pi*random.random_sample(np)
phi = pi*random.random_sample(np)

x, y, z = oovar(size=np), oovar(size=np), oovar(size=np)

S = []
for i in range(1, np):
    tmp = sqrt((x[:i] - x[i])**2 + (y[:i] - y[i])**2 + (z[:i] - z[i])**2)
    tmp.removeAttachedConstraints()
    S.append(sum(1/tmp))
f = sum(S)


x0 = {x : cos(theta)*sin(phi), y : sin(theta)*sin(phi), z : cos(phi)}

p = NLP(f, x0)
p.constraints = [x**2 + y**2 + z**2 == ones(np)]
r = p.solve('ipopt', useSparse=0, options = 'output_file = result.txt, file_print_level = 5')
x_opt, y_opt, z_opt = r(x,y,z)
print(x_opt, y_opt, z_opt)

--------------------------------------------------
solver: ipopt   problem: unnamed    type: NLP   goal: minimum
iter    objFunVal    log10(maxResidual)   
    0  1.493e+03             -15.65
   10  1.071e+03              -0.08
   20  1.056e+03              -2.73
   30  1.055e+03              -3.44
   40  1.055e+03              -3.35
   50  1.055e+03              -4.68
   60  1.055e+03              -4.97
   70  1.055e+03              -6.00
   80  1.055e+03              -9.62
   90  1.055e+03              -8.70
   95  1.055e+03             -11.28
istop: 1000
Solver:   Time Elapsed = 188.04     CPU Time Elapsed = 183.79
objFunValue: 1055.1823 (feasible, MaxResidual = 5.26645e-12)

I guess in Python 3.3. with the LLVM code from Unladen Swallow it will work several times faster.

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

Today I have committed a fix that reduces the time from
Time Elapsed = 188.04     CPU Time Elapsed = 183.79
to
Time Elapsed = 133.63     CPU Time Elapsed = 128.02
(provided you have supplied useSpars=False as I have done in my code)

It could be enhanced even more, but it requires too serious changes in FuncDesigner AD that require toomuch time I can't spend for now.
Also, my NumPy/SciPy in the computer are not linked with AMD ACML (MKL for Intel processors), I guess it could speedup even more.

Regards, D.

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

I've tried running your new vectorized variant with updated OpenOpt and it's great: about 2x faster then previous one (93s). I've also tried building the latest Unladen Swallow and running under it. Now it takes 32s. However I didn't try optimized Numpy (which you've mentioned) too.

Now I understand optimization in Python and vectorization in general much better. Thanks smile

About random numbers: I hope that setting the seed prevents them from being random. That is they should be the same across different script invocations.

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

Try now, with latest SVN changes (r914).
Currently for me it takes
Time Elapsed = 31.33     CPU Time Elapsed = 31.09

Regards, D.

Re: Unexpected inequality constraints in Ipopt with FuncDesigner

Dmitrey wrote:

Try now, with latest SVN changes (r914).

Now I takes 25s in CPython and 29s in Unladen Swallow.