Should I even try to solve this enormous problem with fmincon?

As I am approaching a problem, I realize one way to look at my problem is to consider it as a nonlinear optimization problem. I started looking into the SQP algorithm and fmincon, and I wonder if trying to solve it this way is even a reasonable idea.

I'm trying to minimize f(x), where x is a vector with 4500 elements. f(x) is a double summation of (yij-xj)^2, where i goes from 1 to 50,000 and j goes from 1 to 4500. If it would help, I could use a linear objective function instead, the double summation of (yij-xj).

I then have just over 1 million nonlinear constraints (N*(N-1)/2), where N=1500. These are of the form dmin - sqrt((x1-x4)^2 + (x2-x5)^2 + (x3-x6)^2) < 0. In the solution without constraints, the vast majority of these are satisfied, but I don't know which ones a priori.

Any advice from people who have used fmincon before and understand how its performance scales with problem size / number of constraints?

A brief physical explanation of the problem: I want to find the average (x) over 50,000 configurations the x,y,z coordinates (y) of 1500 particles, with the constraint that the "best fit" configuration cannot have two particles with a separation smaller than dmin.

Answer 1

[Matt J]

You could always try it, though if you do, you should formulate your constraints differentiably

 dmin^2 - ((x1-x4)^2 + (x2-x5)^2 + (x3-x6)^2 )<=0

but I wouldn't expect too much from FMINCON.

It's not so much the large number of constraints, but rather because the feasible set they define looks like it consists of many disjoint subsets that would make getting trapped at local minima very easy. You might try GA in the Global Optimization Toolbox instead.

[Tyler Josephson]

Thanks for the tip. I tried it for two particles, and indeed, the sensitivity of my solution to the initial guess is going to lead to too many problems. For a pair of particles that don't overlap, it gets trapped in a local minimum before it finds the average, unless I use the average as the initial guess. For a pair of particles that overlaps, I cannot use the average as the initial guess, or it will be violating the constraints and trapped in a local minimum.

[Matt J]

Along the lines of Alan's suggestion, maybe you could use one of the Global Optimization Toolbox's multiple start point mechanisms

http://www.mathworks.com/help/gads/how-globalsearch-and-multistart-work.html#bsc9eez

Answer 2

[Alan Weiss]

I believe that you can approach this problem with fmincon if you prepare it carefully. After the reformulation Matt suggested, the gradient and Hessian of both the constraints and objective function are sparse and easily calculable. Use the interior-point algorithm, and sparse Hessian matrices, and fmincon might be able to tackle your problem.

To be clear, you need to set the GradObj, GradConstr, Hessian, and HessFcn options. This example shows a similar setup.

There is no point trying to approach this problem using any other solver, such as GA. If you are afraid of local minima, start the optimization from a variety of initial points.

There are several analyses of disk-packing and sphere-packing problems in the literature, but I do not recall seeing one that used your proposed type of optimization approach. I am not even sure your problem is the same as sphere-packing, but it seems closely related. The Lubachevsky-Stillinger algorithm has been very successful in the past. Unfortunately, my friend Boris Lubachevsky died a few years ago, and I do not know if anyone else has extended his method. Maybe you could give it a try.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation