Global Minima of a Multidimensional Function ??

Hello;

I've recently tried this question in google sci.math.num-analysis DG
with no much luck.
The question appears to be beyond the general interest of the readers
there.

A concise description of the research problem together with my
attempts so far are provided below.

1) The goal is to locate the global minima in a 3D finite domain of a
multidimensional function F(x,y,z).

2) The 3D function F can NOT be expressed analytically or defined in
any manner: expression, equation, formula, or tabulated function, but
its 1st partial derivatives dF/dx, dF/dy, dF/dz can be indirectly
computed at any field point (x,y,z).

3) A 3D uniform orthogonal grid is established. It consists of M*N*K
small 2D panels.
The dimension of the individual elemental panels is approximately of
the same order of magnitude in the 3 directions, and each of the 6
surfaces of 3D elemental "cubical" is small enough to accommodate a
single extremum, if any.

4) Following 2 above, the 3 first partial derivatives are computed at
all grid points. The values of the derivatives are arranged in
increasing order of x, y, and z.

5) I've NOT been able to come up with a numerical procedure for
locating Pm(xm,ym,zm) where the 3 partial derivatives associated with
a minimum (local or global) are diminished simultaneously, i.e.:
dF/dx = 0
dF/dy = 0
dF/dz = 0

6) My approach is that if the 3D function F has a minimum at point P,
then each of the 3 partial derivatives dF/dx, dF/dy, dF/dz has to
simultaneously experience a sign change - to + along the +ve relevant
axis through point P.
{sorry, couldn't type "delta" instead of "d"}

7) I've been able to examine separately each of 6 elemental 2D
surfaces associated with each elemental "cubical", and locate, if any,
a single point (i.e.; vertex) where the applicable 2 of the 3 zero
partial derivatives (item 5 above) are satisfied.
{ The 2D procedure has been extensively checked and validated. }

8) As one would expect, some of the elemental 2D surfaces do not
posses a minimum, while for a given 3D "cubical" potential minimum
could be identified on 1, 2, 3 or even 4 surfaces of its 6 elemental
surfaces.

9) The end result is few vertices V(x,y,z) scattered within the 3D
grid domain:
.....a. About 30 vertices are identified within the 3D orthogonal grid
(consists of 15,000 2D elemental panels):
.....b. some of the vertices are located on the XY elemental panels, 1
vertex per panel, satisfying dF/dx = 0 and dF/dy = 0
.....c. some on the YZ elemental panels, 1 vertex per panel,
satisfying dF/dy = 0 and dF/dz = 0
.....d. some are on the 2D XZ elemental panels, 1 vertex per panel,
satisfying dF/dx = 0 and dF/dz = 0
.....e. Only few of these vertices are likely to be associated with the
global minima, and the rest are most likely locals.
.....f. If it helps, I'd be glad to provide the coordinates (x,y,z) of
the 30 or so vertices and their property values (item 10 below) }

10) With respect to the "global" minima, fortunately, I can evaluate a
"property" at each located vertex, the value of which is inversely
proportional to the value of the unknown function F at the same
vertex.
.....g. In other words, the greater the "property" value is, the closer
the vertex to the global minima.
.....h. The 3D function F has at least one minimum in the domain, and F

11) The question is:
How to mathematically/numerically manipulate those vertices (item 9
above) to locate a single point Pm(xm,ym,zm) {global minima, with
reasonable figure of merit) where the 3 partial derivatives are
satisfied simultaneously.

Your help would be greatly appreciated.

Thank you kindly.
Monir

 

CORRECTION:
===========

Sorry! Item 9.a of my post should read:

9) The end result is .....
.....a. Total 87 vertices are identified within the 3D orthogonal grid
(of 15,000 2D elemental panels), each satisfying 2 of the 3 conditions
for minimum :
54 vertices on y-z elemental surfaces, where dF/dy=0 & dF/dz=0
16 on x-z elemental surfaces, where dF/dx=0 & dF/dz=0
16 on x-y elemental surfaces, where dF/dx=0 & dF/dy=0

My apologies!
Monir

 

I wonder why you didn't try Roger Kearfotts INTBIS which is available
from http://www.netlib.org/toms.
(computes all zeroes of a function on a box)
If I understood your question right, when you try to find the global
minimizer of a real function of three variables on a box by identifying
all zeroes of the gradient. you can evaluate the gradient at given points,
and then evaluating an optimality criterion related somewhat to your
objective function which cannot be evaluated directly. also you state
that your grid boxes "are small enough to accomodate a single extremum"
does that really mean that you know upper bounds for the inverse of the
Hessian?
maybe a look on Arnold Neumaiers global optimization page which also
covers the zero finding problem is also helpful, if pointwise function
values are all you can provide.
http://www.mat.univie.ac.at/~neum/globt.html
hth
peter

 

Hi Peter;

Will review netlib INTBIS procedure /toms/681.
Please keep in mind that the 3D function F is unknown (equation &
values), and although values of grad(F) is available at ALL grid
points, it can NOT be expressed analytically or defined in any manner.

Regards.
Monir

 

(2nd posting of my earlier reply of Sep 18, 2009)
Hi Peter;

You're absolutely correct!
1) The package INTBIS (/toms/681) contains all what is needed to use
it. I didn't appreciate earlier that it had been developed some time
ago, modified, added to, etc., and unavoidably some old comments had
been left in.
So when its documentation states, for example: "INTBIS includes the
test input file TOMS86SH.DT1, the test configuration file CONFIG.BIS
and the corresponding output file TOMS86SH.OU1", I looked for those
specific files, and didn't take it to mean that the example input/
output files are included at the end of the source code with no name
designation!

2) Let me first describe my interpretation of your suggestion.
Basically, we will treat the gradient vector (3 comp) known at all the
3D grid points as 3 (independent) functions, and use INTBIS to find
the points within the grid where the function vector (dCpX,
dCpR,dCpTH) vanishes.
This netlib procedure, if successful, would locate the points pi(xi,
ri, thi) of zero function vector, i.e.; real roots (if any), where the
unknown function P takes on a maximum or a minimum value.
Among those extrema, one would need (later) to identify the global
minima point pm(xm,rm,thm) of the function P.
Will address that when we get there!!

3) I've got the INTBIS algorithm working properly with the provided
test cases (see item 4 below). The output correlates extremely well
with the results included in package.
Some exceptions, however, regarding the residual vector at some
roots. My run appears to produce more realistic residual results of
the order 0.10D-09 compared with the documentation-provided residual
of the order -0.19D-01 at the same root (box 2 root, Test problem no.
14.).
(See item 5 below)

4) To achieve that, I had to make some changes to the code, including
modifications to few variable array declarations, populating input
data file FORT.20 (connecting the read device UNITC), activating the
2nd m/c constants set in D1MACH(), and activating the 1st set in I1MACH
().
My selection of the above m/c constants sets was primarily based on my
guess of what the *most applicable* sets are.

5) Surprisingly, I haven't made any changes to SIMINT() or to ERRHND
().
The latter routine has no executable statements. Its Call ERRSET()
statement is commented out, and the routine ERRSET() is not included
in the package!
(This might be contributing somehow to my observation in item 3.
above.)

6) My next task was to re-examine Sub POLFUN(). It is of no use to me
as you kindly suggested.

SUBROUTINE POLFUN(N,X,SX,FVAL,SCALF)
DOUBLE PRECISION X(2,N), SX(N)
DOUBLE PRECISION FVAL(2,N),SCALF(N)

C N is the number of equations and variables.
C (INPUT)

C X X(1,I) is the left endpoint of the I-th coordinate interval,
C and X(2,I) is the right endpoint of the I-th coordinate
C interval of the box, for I between 1 and N.
C (INPUT)

C SX contains a point N-vector at which a point function value is
C to be computed.
C (INPUT)

C FVAL Upon return, ( FVAL(1,I), FVAL(2,I) ) will be an interval
C which contains the range of the function over the box X.
C (OUTPUT)

C SCALF Upon return, SCALF will contain an approximate vector
function
C value at the point SX.
C (OUTPUT)

7) I've written (and successfully tested) a 3D spline interpolation
double precision routine POLIN3d(), which returns the interpolated
function value at any point (X1,X2,X3) within the 3D initial grid.

SUBROUTINE POLIN3d(x1a, x2a, x3a, YA, kx, mr, nth,
1 x1, x2, x3,
2 Z )
C***********************************************************************
! Given arrays x1a (length kx), x2a (length mr) and x3a (length nth)
of independent
! variables, and a kx by mr by nth array function values YA ,
tabulated at the grid
! points defined by x1a, x2a and x3a; and given values x1, x2 and x3
of the
! independent variables; this routine returns an interpolated function
value Z.

8) Back to INTBIS Sub POLFUN() item 6. above.
My approach is to keep the routine name and replace its contents with
Call POLIN3d(); 9 calls in total. Something like:

SUBROUTINE POLFUN(N, X, SX, FVAL, SCALF)
C***********************************************************************
! This modified routine is exclusively for N = 3, and DOES NOT use
! or reference to poly coeff, deg, terms, etc.
!.........................................................................................
double precision X(2,N), SX(N)
double precision FVAL(2,N),SCALF(N)
!.........................................................................................
! 3 calls to evaluate functions at SX(1:N), N=3
Call POLIN3d(x1a, x2a, x3a, dCpX, kx, mr, nth,
1 SX(1), SX(2),SX(3),
2 SCALF(1) )

Call POLIN3d(x1a, x2a, x3a, dCpR, kx, mr, nth,
1 SX(1), SX(2),SX(3),
2 SCALF(2) )

Call POLIN3d(x1a, x2a, x3a, dCpTH, kx, mr, nth,
1 SX(1), SX(2),SX(3),
2 SCALF(3) )
...........................................................................................
! 3 calls to evaluate functions at LEFT endpoint of interval X
(1,1:N), N=3
Call POLIN3d(x1a, x2a, x3a, dCpX, kx, mr, nth,
1 X(1,1), X(1,2),X(1,3),
2 FVAL(1,1) )

Call POLIN3d(x1a, x2a, x3a, dCpR, kx, mr, nth,
1 X(1,1), X(1,2),X(1,3),
2 FVAL(1,2) )

Call POLIN3d(x1a, x2a, x3a, dCpTH, kx, mr, nth,
1 X(1,1), X(1,2),X(1,3),
2 FVAL(1,3) )
...........................................................................................
! 3 calls to evaluate functions at RIGHT endpoint of interval X
(2,1:N), N=3
Call POLIN3d(x1a, x2a, x3a, dCpX, kx, mr, nth,
1 X(2,1), X(2,2),X(2,3),
2 FVAL(2,1) )

Call POLIN3d(x1a, x2a, x3a, dCpR, kx, mr, nth,
1 X(2,1), X(2,2),X(2,3),
2 FVAL(2,2) )

Call POLIN3d(x1a, x2a, x3a, dCpTH, kx, mr, nth,
1 X(2,1), X(2,2),X(2,3),
2 FVAL(2,3) )

9) The major problem now is what to do with the read poly no. of
terms, coeffs, powers, etc.
The algorithm uses this data throughout! ... e.g.; in routines Sub
POLJAC(), Sub POLJSC(), etc. And if not, it passes the data in common
blocks to other routines.
It's a major headache!!

10) Changing the I/O of the procedure would not only be a major
undertaking but also a dangerous proposition!
A quick fix: How about having the UNITI input data to include a single
and bogus poly non-linear equation with ALL coefficients specified as
0.0 for all terms ??
I don't believe it will work. Will keep trying, and will keep you
posted if interested!

Would very much appreciate any comments or suggestions you might have.

Regards.
Monir

PS. Please let me know if you prefer that we correspond by emails so
not to take too much space here. Once a successful resolution to the
issue is reached, I'd be glad to post here (sci.math.research) a
concluding summary of the solution.