Ask Your Question
2

Black box numerical optimization of a convex function

asked 2014-08-19 05:43:57 -0500

Thomas gravatar image

updated 2014-08-19 09:23:47 -0500

I would like to numerically optimize a convex function without using any of its internals, that is, a python function. Is there an implementation of this in sage? I tried sage.numerical.optimize.find_local_minimum but it fails internally with a

TypeError: unable to find a common ring for all elements

My function to optimize is in fact the second largest eigenvalue of a matrix that is constructed from one parameter. Here it is in all gory detail:

def binaryStrings(k):
    for i in range(2**k):
        s = str(bin(i))[2:]
        yield "0"^(k-len(s)) + s;

def hammingDistanceInt (a,b):
    return bin(a^^b).count("1")

def cubeAdjacency (k):
    # Note: We would like to use QQ instead of CDF, but this will
    # result in inexact computations failing due to rounding later
    # (Yes it does!)
    # https://groups.google.com/forum/#!topic/sage-support/3IIbu0OBJX4
    return matrix(CDF, 2^k, 2^k, lambda i, j: 1 if hammingDistanceInt(i,j)==1 else 0);

def hemmeckeFiberAdjacency (k, p1):
    p2 = (1 - 4*k*p1)/2;
    M = p1 * block_diagonal_matrix (cubeAdjacency(k),cubeAdjacency(k));
    M[0,2^k] = p2
    M[2^k,0] = p2
    rowsums = M * vector ([1]*2^(k+1))
    M += diagonal_matrix(vector([1]*2^(k+1))-rowsums);
    return M

def slem(k,p):
    e = map(abs, hemmeckeFiberAdjacency(k,p).eigenvalues())
    e.sort();
    return e[-2];

def slem4(p): return slem(4,p)
plot(slem4, [0,1/16]) # Works fine
sage.numerical.optimize.find_local_minimum(slem4, [0,1/16]) # Type error
edit retag flag offensive close merge delete

Comments

1

In order to understand where the error comes from, could you please provide the construction of the function to optimise ?

tmonteil gravatar imagetmonteil ( 2014-08-19 07:48:58 -0500 )edit

1 answer

Sort by ยป oldest newest most voted
3

answered 2014-08-22 19:34:41 -0500

tmonteil gravatar image

Here is the problem: if you look at the Traceback, you can see that the problem is in hemmeckeFiberAdjacency, when you compute the diagonal_matrix. You can also see that this function tries to find a common ring for the elements you give them.

So let us add some bugtracking in your code by replacing:

M += diagonal_matrix(vector([1]*2^(k+1))-rowsums);

by:

try:
    M += diagonal_matrix(vector([1]*2^(k+1))-rowsums);
except Exception:
    print "bug", p1, p1.parent()

so that when the error will appear, we will be able to see which p1 caused the problem.

Indeed, we get another error : AttributeError: 'numpy.float64' object has no attribute 'parent'. Here is the explanation : when you call sage.numerical.optimize.find_local_minimum, Sage uses scipy.optimize which uses numpy.float64 floating point numbers. Those numbers have no parent (a ring to live in) in Sage.

To fix your problem, you can just ensure that the p1 that is given to hemmeckeFiberAdjacency can be correctly handled by Sage. It suffice to transform it into an element of the Real Double Field (which are Sage floating point numbers with the same precision as numpy.float64 elements) by adding the following line at the beginning of your function hemmeckeFiberAdjacency.

p1 = RDF(p1)

Now you will get:

sage: sage.numerical.optimize.find_local_minimum(slem4, 0,1/16)
(0.993614695234, 0.042890047690081409)

Which corresponds to what you saw on the picture.

edit flag offensive delete link more

Comments

Thanks a lot for the detailed explanation!

Thomas gravatar imageThomas ( 2014-08-31 10:47:03 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2014-08-19 05:43:57 -0500

Seen: 94 times

Last updated: Aug 22 '14