R1NNOOY KAN Econometric Institute, Erasmus Llniver~itv Rotterdam, The Netherlands

24  Download (0)

Full text

(1)

Mathematical Programming 37 (1987) 184 207 North-Holland

H I T - A N D - R U N A L G O R I T H M S F O R T H E I D E N T I F I C A T I O N O F N O N R E D U N D A N T L I N E A R I N E Q U A L I T I E S

H.C.P. BERBEE

Department ~1 Mathematical Statistics, Centre./~r Mathematics and Computer Science, Amsterdam,

"Hie Netherlands

C.G.E. B O E N D E R

Econometric Institute, Erasmus Universi O, Rotterdam, "Hie Netherlands

A.H.G. R1NNOOY KAN

Econometric Institute, Erasmus Llniver~itv Rotterdam, The Netherlands

C.L. S C H E F F E R

Department ol Mathematics and In lT~rnlatics, De!t7 University of Technology, V/he Netherlands

R.L. SMITH

Department r~/'lndustrial and Operations Engineering, University of Michigan, Ann Arbor, Michigan, USA

J. T E L G E N

Van Dien + ('o., Utrecht, The Netherlands

Received 11 February 1985 Revised I1 June 1986

Two probabilistic hit-and-run algorithms are presented to detect nonredundant constraints in a full dimensional system of linear inequalities. The algorithms proceed by generating a random sequence of interior points whose limiting distribution is uniform, and by searching for a nonredundant constraint in the direction of a random vector from each point in the sequence. In the hypersphere directions algorithm tile direction vector is drawn from a uniform distribution on a hypersphere. In tile computalionalb superior coordinate directions algorithm a search is carried out along one of the coordinate vectors. The algorithms are terminated through the use of a Bayesian stopping rule. Computational experience with the algorithms and the stopping rule will be reported.

Key words: System of linear inequalities, redundancy, probabilistic hit-and-run algorithms, uniform interior points, Bayesian stopping rule.

1. Introduction

The problem o f recognizing redundant linear inequalities (i.e., inequalities that can be deleted from a system without changing its set o f feasible solutions) is o f obvious computational importance. Methods to eliminate such constraints have been proposed by many researchers. Most o f them are based on the simplex method (e.g. Thompson et al., 1966; Lisy, 1971; Gal, 1975; Telgen, 1979); sometimes, the

184

(2)

H.C.P. Berbee et al./ Hit-and-run algorithms to ident~y]' nonredundant c~mstraints 185

simplex method is in fact invoked to e n u m e r a t e a l l the extreme points of the feasible solution set (Balinski, 1961; Mattheis, 1973). A detailed computational comparison (Karwan et al., 1983) reveals that all these approaches suffer under the administrative burden of maintaining a complete updated simplex tableau. As a result, they require a computational effort that is hardly compensated for by the subsequent speed-up of calculations carried out on the reduced set of inequalities.

An alternative for these time consuming procedures is provided by

heuristic

methods. Acting as

preprocessors

on the original problem data (Brearly et al., 1975;

Bradley et al., 1980), many of these heuristics have found their way into commercial mathematical p r o g r a m m i n g packages, for example under the name of the

REDUCE

option, Of course, the big drawback of these fast procedures is that there is no guarantee that all redundancy present in the system will be identified.

An attractive compromise between the two above approaches is provided by the

probabilistic preprocessors

(cf. Rabin, 1976) that are studied in this paper. The aim of these preprocessors is to identify n o n r e d u n d a n t constraints rather than redundant ones. Starting from some interior point (assumed given) of the polyhedron defined by the system (assumed to be b o u n d e d and full dimensional), the preprocessor generates a random sequence of interior points: in each successive point, a search is carried out in a r a n d o m direction and the constraints first encountered in that direction and its negation are identified as being nonredundant (cf. Theorem 1). In the

hypersphere directions method,

the random direction is generated from a uniform distribution on a hypersphere (Boneh and Golan, 1979; Smith, 1980; Boneh, 1983);

in the

coordinate directions method,

it is chosen with equal probability from the coordinate direction vectors and their negations (Telgen, 1980). In both cases, the r/ext interior point is generated randomly from a uniform distribution over the line segment connecting the two points where the direction vector and its negation intersect with the polytope. Methods of this type, for obvious reasons, are also referred to as

hit-and-run algorithms.

(In Smith (1984) the term

symmetric mixing algorithms

is used.)

For the hypersphere directions method, it is shown in Smith (1984) that the sequence of interior points has a

limiting distribution

which is uniform over the interior of the polytope. We give a new proof for this result and extend it to the (much more difficult) case of the coordinate directions method in Section 3. These results imply that both methods are

asymptotically correct

in the sense that each nonredundant constraint will be identified with probability one as the number of iterations increases. Indeed, the computational experiments in Karwan et al. (1983) indicate that the great speed at which interior points can be generated turns such a method into a very attractive preprocessing device of great practical value,

This good practical performance provides part of our motivation to consider the important question of when to terminate a hit-and-run procedure; under the assump- tion that the true n u m b e r of nonredundant constraints is unknown, the answer is not obvious. To cope with this problem we will prove in Section 4 that the theorems of Section 3 imply that asymptotically each

facet

of the boundary of the feasible

(3)

186 H.C.P. Berbee et aL / Hit-and-run algorithms to idenHl~v nonredundant constraints

region (or, equivalently, each n o n r e d u n d a n t constraint) has a fixed probability o f being hit. In Boneh a n d G o l a n (1979) this result is applied to determine the expected n u m b e r o f iterations that the hypersphere directions algorithm will need to identify all n o n r e d u n d a n t constraints, u n d e r the additional a s s u m p t i o n s that (i) the success- ive interior points are statistically i n d e p e n d e n t , (ii) the a s y m p t o t i c hitting prob- abilities are all equal, a n d (iii) each constraint is n o n r e d u n d a n t . In practice, however, the hitting probabilities are very different: in systems with less than 100 constraints in 20-dimensional space, for example, we frequently observed empirical probabilities smaller than 10 4. Boneh and G o l a n regard their result as a I o w e r b o u n d for the n u m b e r o f iterations, and provide an u p p e r b o u n d u n d e r the a s s u m p t i o n that the (ordered) probabilities decrease arithmetically. H o w e v e r , under the a s s u m p t i o n (i) we will show in Section 4 that the generalized multinomial model studied in Boender (1984, C h a p t e r 3) applies, so that we can obtain a more accurate Bayesian stopping criterion for both hit-and-run algorithms, in which the user m a y i n c o r p o r a t e a priori information on the equality of the hitting probabilities and on the true n u m b e r o f n o n r e d u n d a n t constraints, instead o f a s s u m i n g (ii) and (iii). In particular, the hit-and-run algorithms can be terminated if the optimal Bayesian estimate o f the n u m b e r o f yet u n d i s c o v e r e d n o n r e d u n d a n t constraints is 0.

O u r experimental results are presented in Section 5; c o n c l u d i n g remarks and possible extensions to optimization p r o b l e m s (e.g. linear p r o g r a m m i n g ) are the subject o f Section 6.

2. The hit-and-run algorithms

C o n s i d e r a feasible region S defined by a system o f linear inequalities

a [ x < ~ b i ( i = 1 . . . . , m ) (1)

with x ~ ~a and ]] a, ]] = 1 (i -- 1 , . . . , m ). We will assume that S is b o u n d e d , n o n e m p t y and o f full dimension, so that S is a polytope that contains interior points for which the inequalities (1) are all satisfied as strict inequalities. A redundant constraint is defined as an inequality which may be d r o p p e d from the system (1) without changing the feasible region S. A Jacet is defined as the interaction o f a n o n r e d u n d a n t constraint with the b o u n d a r y o f S o f d i m e n s i o n n - 1.

O u r probabilistic hit and run preprocessors to investigate r e d u n d a n c i e s within system (1) are based on a search from an interior point X in the direction of a vector v with tlvt[= 1. Let us d e n o t e the straight line passing t h r o u g h X in the direction v by X + A v (AcIR). Then it is immediate that the value o f A at the intersection point with the i-th h y p e r p l a n e aV;x = b, is equal to Ai = (bi - a V ; X ) / a r v . A c c o r d i n g to the following theorem, the constraints hit first in the positive (A > 0) a n d negative (3, < O) direction can be declared to be n o n r e d u n d a n t .

Theorem 1. I f

r a a rgm}n{,L ]3,, > O} (2)

(4)

H.C.P. Berbee et a l . / Hit-and-run algorithms to ident(/.i' nonredundant constraints 187

and

s ~ argmax{A~ I A, < 0} (3)

are unique, then the constraints a T x <~ br and a l,-x <~ b, are nonredundant.

Proof. We o n l y offer a p r o o f for the case c o r r e s p o n d i n g to (2). It suffices to s h o w (Telgen, 1979) that u n d e r the c o n d i t i o n s o f the t h e o r e m , t h e r e exists a p o i n t X "

such t h a t

aV,.X"> b,., (4)

a T X " ~ b i ( i = 1 . . . m , i # r ) . (5)

O b v i o u s l y A , = ( h i - - a T X ) / a T v > O for an i n t e r i o r p o i n t X if[ a T v > O . H e n c e , for X ' = X + A,.v we have, f o r i # t;

T , b , . - T a r X T bi - a T X T

a, X = a T X ~ r a, v < a / ' X -~ - ai v = b,. (6)

a r t 2 U/ L~

Since a T X '= b~, there exists an e > 0 such that X .... X ' + ev satisfies (4) a n d (5).

If (2) a n d (3) are n o t u n i q u e , then e i t h e r i d e n t i c a l c o n s t r a i n t s h a v e b e e n hit, or a n i n t e r s e c t i o n o f c o n s t r a i n t s has b e e n hit. The f o r m e r p o s s i b i l i t y is a s s u m e d not to o c c u r from now on. T h e l a t t e r p o s s i b i l i t y can o n l y o c c u r with p r o b a b i l i t y 0, a n d t h e r e f o r e d i s r e g a r d e d .

The a l g o r i t h m s to be p r e s e n t e d e x p l o i t the a b o v e t h e o r e m by s e a r c h i n g in the d i r e c t i o n o f a r a n d o m v e c t o r v" from each p o i n t X " o f a r a n d o m sequence o f i n t e r i o r points. F o r the hypersphere directions algorithm ( B o n e h a n d G o l a n , 1979; Smith, 1980; B o n e h , 1983], a d i r e c t i o n v e c t o r is d r a w n in each p o i n t X" o f the s e q u e n c e from a u n i f o r m d i s t r i b u t i o n on a unit h y p e r s p h e r e with c e n t r e X". The coordinate directions method (Telgen, 1980) g e n e r a t e s o n e o f the unit c o o r d i n a t e vectors or t h e i r n e g a t i o n as d i r e c t i o n vector. Both a l g o r i t h m s c h o o s e the ( n + l ) - t h i n t e r i o r p o i n t un(/brmly on the l i n e s e g m e n t c o n n e c t i n g the p r e v i o u s two hitpoints, i.e. the p o i n t s w h e r e the line X " + Av" intersects the b o u n d a r y o f S. The c h o i c e o f a s t o p p i n g c r i t e r i o n is d i s c u s s e d in Section 4.

The a l g o r i t h m s c o n s i s t o f the f o l l o w i n g steps.

Step 0: F i n d an i n t e r i o r point X ~ Set n := 0.

Step I: Hypersphere directions algorithm: G e n e r a t e a d i r e c t i o n v e c t o r v" from a u n i f o r m d i s t r i b u t i o n o n a unit h y p e r s p h e r e with c e n t r e X " . Coordinate directions algorithm: G e n e r a t e a d i r e c t i o n v e c t o r v" with equal p r o b a b i l i t y from one o f the d c o o r d i n a t e vectors a n d t h e i r n e g a t i o n s .

Step 2: D e t e r m i n e bi -- a ~ X "

A,: T , ( i = 1 . . . m), (7)

a i u

(5)

188 H.C.P. Berbee et al. / Hit-and.run algorithms to identi[i" nonredundant constraints

A +:= rain {A~IA~>0}, (8)

A := max {A, I A , < 0 }. (9)

Declare the constraints corresponding to the indices for which the minimum (8) and the maximum (91) are attained to be nonredundant.

Step 3: Generate u from a uniform distribution on [0, 1] and set

X " + ~ = X " + ( A - + u ( A ' - A ))v". (10)

Step 4: Set n := n + 1 and go to Step 1, unless a stopping criterion is satisfied.

We conclude this section with a comparison of the computational efficiency of the two algorithms. It is easily verified that the determination of the intersection points of a direction vector with the m hyperplanes (cf. (7)) is the most time consuming part of both algorithms. Since S c ~J, the computation of one intersection point for the hypersphere directions algorithm requires O(d) time, which implies that it requires O(rod) time to evaluate the intersection points of a straight line with all the m hyperplanes. The major advantage of the coordinate directions algorithm is that since only a single coordinate is changed when moving from one interior point to the next, no more than 2 multiplications are needed to update A~

( i = 1 , . . . , m) (cf. (7)). Hence, for the coordinate directions algorithm the computa- tion of the intersection points requires O(m) time, rather than O ( m d ) .

3. The uniform limiting distribution of the interior points

In this section we will prove that for both hit-and-run algorithms the random sequence {X'}i" o of interior points converges to the uniform distribution /~ on S, for each possible starting point in its interior S ~'. These results are an essential part of the justification of our application of the generalized multinomial model, which underlies the stopping criterion proposed in Section 4. We remarked earlier that for the hypersphere directions algorithm this result has already been proved in (Smith, 1984). However, our proof for this case is new, and it provides an introduction to the p r o o f for the coordinate directions algorithm.

For the remainder of this section we fix a starting point x ~ ~ and a Borel set B ~ S with I*(B")> 0, and without loss of generality we assume that the Lebesgue measure of S is equal to 1.

From the description of the algorithms in the previous section it is immediate that for all n~[~ ~,and x l , . . . , x " c S

P r { X " + ' c B ~ I X ~ x" . . . X" = x"} = P r { X " ~ ' c B~ X'' = x " } , (11 and that

Pr{X"* 1 ~ B " I X " = x ~'} -- Pr{X' ~ B" I X " = .x"}. (12 Hence, for both algorithms the sequence of interior points defines a M a r k o v chain with a stationary transition probability.Junction and continuous state space S.

(6)

H.C.P. Berbee et a l . / Hit-and-run algorithms to identil)' nonredundant constraints 18 9

Theorem 2. The hypersphere directions algorithm generates a sequence of interior points whose limiting distribution is unffbrm on S, i.e.:

lira Pr{X" c B"] X ~ x ~ = tz(B~ (13)

Proof. A c c o r d i n g to T h e o r e m 7.1 in Orey (1971) it is sufficient to prove P r o p o s i n o n s (i) and (ii).

(i) tx is invariant, i.e. if a current interior point is uniformly distributed, then the next interior point is uniformly distributed as well:

, Pr{X' 6 B ~ ~ x} dx = ~ ( B ~ (14)

(ii) The M a r k o v chain is ~-recurrent, i.e. B ~ is visited at least once with probability 1"

Pr{3i6[~+: X ' ~ B ~ 1 7 6 1 7 6 1. (15)

To prove (i) and (ii), we first derive the expression for the transition densio, Jhnction p : SO• S ~ R + w {0}. For x, y c S ~ (x # y), let Hy d e n o t e a d - d i m e n s i o n a l h y p e r c u b e with centre y and volume h a, oriented along the ray from x to y. Then the transition density o f x, y is defined as

p ( x , y ) = l i m P r { X ' e H,,IX ~ x}

~'~" h a (16)

D e n o t e the surface area o f a d - d i m e n s i o n a l h y p e r s p h e r e with radius r by c(r), and let m ( x , y ) be the diameter o f S m e a s u r e d along the ray from x to y. Then, since the h y p e r s p h e r e directions algorithm in x generates its direction vector uniformly on the (unit) h y p e r s p h e r e with centre x, a n d chooses the next interior point uniformly on the intersection o f this direction vector with S we have (cf. Figure 1):

P r { X ' c H , t X ~ 1 2h ~ ' h 2

lira = lira ha - . ( 1 7 )

,,- c({ly-xll) m(x,y)

c({ly-x{I)m(x,v)

h 1 o h d

Hence,

fsPr{Xl~B~176 I p(x,y)dydx

= ,3" c ( I J Y - x l l ) - rn(x, y) dy dx

u" s c ( I ] x - y l l ) 9 re(y, x) dx dy

= f Pr~X' ~SIx~ f

d y = . ( B ~

Bo B t~

(18) which proves (i). We observe that for this p r o o f it is sufficient that a transition density function exists, and that it is s y m m e t r i c in its arguments.

(7)

190 H.C.P. Berbee et at./ Hit-and-run algorithms to idenlit ), nonredundant constraints

Fig. 1.

To prove (ii), let

&

rs = sup [ ] y - x [ [ .

.x,v~ 5

(19)

Since S is assumed to be b o u n d e d , rs is finite. Hence

P r { X ' ~ B ~ " = x " } = ~ p ( x " , y ) d y 3 B"

I 2 f 2

x I I ) m ( x , d v ~ - - d r

~,, c / [ l y - " " ~' y) " t~" c(rs)rs 2/~(B ~

- (20)

c( "s ) rs

Thus, since the lower h o u n d (20) is independent o f x ~ we have

Pr{B i e ~+: X i e B ~ ~ = x ~ = 1 - Pr{X ~ ~ B ~ V i e [~+lX ~ = x ~

= 1 - l i r a P r { X ~ B ~ = 1 , . . . , j ] X " = x ~}

j ~ : x ,

In a d d i t i o n to x ~ and B ~ we define H ~ to be a h y p e r c u b e with p . ( H ~ which is fully c o n t a i n e d in & a n d whose edges h, ( i = 1 . . . d) are oriented along the c o o r d i n a t e axes e~ (i = 1 . . . . , d).

Theorem 3, The coordinate directions algorithm generates a sequence o f interior points whose limiting distribution is uniform on S:

lira P r { X " c B ~ ~ x ' ) : . = / , ( B~ (22)

n ~ , x

Proof. A n a l o g o u s l y to the p r o o f o f T h e o r e m 2 we proceed by proving Propositions (i) a n d (ii).

We note, however, that given the location x o f the n-th interior point, X "§ will be c o n t a i n e d in the set o f c o o r d i n a t e axes t h r o u g h x. This set is o f g - p r o b a b i l i t y 0,

(8)

H . C . P . B e r b e e e t a l , / H i t - a n d - r u n a l g o r i t h m s t o i d e n t ( l ; v n o n r e d u n d a n t c o n s t r a i n t s 1 9 1

which implies that for the coordinate directions algorithm the l-step transition density function is not defined. Hence, to prove (i) we cannot apply the simple method of proof of Theorem 2.

To prove Proposition (i) it is sufficient to show that (cf. (14))

fs Pr{X~ e H~176

x} dx = # (H~ (23)

Define T~ c S as step if a move is

the set of points starting from which H ~ can be reached in one made along e~ ( i = 1 .. . . . d). Then (cf. Figure 2):

t P r { X ' c H ~

.g

= ~" I P r { X l c H ~ Pr{X'cH~176

i = I T, H I' H u

= ~ I Pr{X'eH~177176

i = I T, t~l ~

't I

+ v Pr{X ~

i = 1 , H u

c H ~ v ~ •

~ x} dx

= ,., | Pr{X'~H",v~177 lX"=x}dx.

f (24)

i = 1 d l ;

Choose y S T~, and define FI. as the intersection o f S with the straight line through y along direction e,. Then, for all x c F',.,

Pr{X'c H ~ v" • l x ~ Iit,,11 (25)

dllF;ll'

which does not depend on x. Hence,

f Pr{X ~ r ~ !) ~ = •

~

f

- -

v; ./ v~.

if xi denotes the i-th component of x, then tlh'lJ dx,

Ijh, ll ~ jjh, ij

a

FI.I -allF;ll IICII= d

(26)

F i g . 2 .

(9)

1 9 2 H.('.P. Berbee et al./ Hit-and-run algt~rithm~" to ident([y nonredundant constraints

Thus, integration over all remaining components j # i gives

f P r { X ' c H

~176177

E;L, IIh~ll (27)

T, d

and, substitution of (27) in (24) yields

PF{XI c

n~

dx = 11 IIh, II = /x(H"), (28)

which proves (i).

From (21) it should be clear that the p r o o f of Proposition (ii) follows easily if

inf Pr{.If I c B~ ~ x} > 0. (29)

However, as only moves along the direction of the coordinate vectors are allowed, (29) is obviously not true. Furthermore, one can easily devise polytopes S for which (29) does not hold if instead of one perpendicular step any finite n u m b e r n of such steps may be taken, so that this problem cannot be remedied by considering n step transition probabilities either.

A diflerent approach is therefore necessary to prove Proposition (ii). Our p r o o f will be based on Proposition 5.1 in Orey .(1971), in which the initial starting point is assumed to be a random variable with distribution rt rather than being fixed at

x".

Let Pr,,{C} stand for the probability' of an event C, given that X ~ is distributed according to rl, and choose a real scalar c2'> 0. Then Orey's Proposition states that, if a set A c S exists (possibly depending on B" and c2 ~) such that

and

then

i ~- o

Pr,{{X }, :o~ A infinitely' o f t e n } > 1 - c~ (30)

P r , ~ { ( X ' } i ' o c B " i n f i n i t e l y o f t e n } > l c2 ~. (32) This implies in turn that

Pr,.,{{X'}~oc. B ~ infinitely often} = 1, (33)

and a fortiori

Pr,~{3iE ~ : X ' e B~ = 1. (34)

Let the set A~ be defined as the subset of those points in S whose distance to the b o u n d a r y of S is greater than e. > 0. Then, with Orey's Proposition in mind, we will

inf P r { 3 i ~ N * : X ~ c B ~ (31)

(10)

H.C.P. Berbee et aL / Hit-and-run algorithms to ident![~v nonredundant constraints 193

prove the following statements l and II, u n d e r the additional a s s u m p t i o n that r/ is b o u n d e d from above b y / z , i.e. there exists a constant c > 0 such that r / ( E ) <~ c/z ( E ) for all Borel sets E c o n t a i n e d in S.

I. There exists an e ( B ~ such that for all e ~< e(B~ A = At- satisfies (31).

II. There exists an e ( a ~ such that for all s ~< e(a~ A = A~ satisfies (30).

If I and II are true, then A = At with e -- min{~(B"), s ( a ~ simultaneously satisfies (30) and (31), which then proves (34). We then have proved (ii) for the case that the starting point follows an initial distribution which is b o u n d e d from above by /~. It then still remains necessary to extend (34) to the case that the process is started in x ~

We will first make that final step. The p r o o f consists o f first p e r f o r m i n g n steps o f the procedure, starting from x ~ The distribution o f X" then will serve as the initial distribution r/, and the result desired is o b t a i n e d by applying (33), assuming that I and I1 are satisfied.

By the Lebesgue D e c o m p o s i t i o n T h e o r e m (Ash, 1972) the n step transition probability distribution can be d e c o m p o s e d into a singular part, o,, say, whose probability mass is c o n c e n t r a t e d on a set with /z-probability 0, and an absolutely continuous part, w,, say, with an integrable function p"(x~ . ) so that for all Borel sets E in S

,o,,(E) = r p ( x , v) dr. " " . (351)

J E

Clearly the M a r k o v chain is c o n c e n t r a t e d on a set o f 0 /z-probability if not all d distinct c o o r d i n a t e vectors have already been generated as r a n d o m direction. If n ~>d it is easily verified that the probability o f this event does not exceed d ( 1 - ( 1 / d ) ) " , so that also

v,,(S)<~d 1 - (36)

for n ~> d. N o w define

a3,,& w, (37)

1 -v,,(S)"

Then a3,,(S)--1, so that o3,, is a probability distribution with a probability density function. F r o m the description of the c o o r d i n a t e directions algorithm and from (35) it follows that oS,,(E) is determined by n successive 1-dimensional integrations o f constant functions over line segments t h r o u g h x ~ ' , X ~ , . . . , X " ~ whose function value is the inverse o f the length o f the line segments. Given any n and o u r starting point x", these line segments are b o u n d e d from below by a positive constant, so that o3,, is b o u n d e d f r o m above b y / z . Hence, assuming 1 and I1, (33) is satisfied for

rt = oOn. Furthermore, for each Borel set E c o n t a i n e d in S

(11)

194 H.C.P. Berbee et al./ Hit-and-run algorithms to ident(ly nonredundant constraints

1 - o,,(S)

Pr{X" c E IX ~ x ~

= w,,(E) + v,,(E) = w,,(E) ~- v,,(E)

l - o , , ( S )

> o~,,(E)(1 - ~,,,(S)). (38)

Hence, by (36) and (33) we have, for all n ~ d, Pr{{X~};~_o ~ B ~ infinitely often IX ~ = x ~

= [ Pr{X" ~ dy IX" = x ' ) Pr{{X~}L,, c B ~ infinitely often

I X ~ = y}

3 S

> ( 1 - o , , ( S ) )

[

tb,,(dy) Pr{{X*}Z,,< B ~ infinitely often ] X ~

.1,

= ( 1 -

o,,(S)) Pr,a,,{{X'}~_oC B ~

infinitely often}

~> Pr,a,,{{X };=~,c B ~ infinitely often}

( -•

= l - d 1 d ! , (39)

so that

P r { { X ' } ~ , , c B ~ infinitely o f t e n l X " = x u} = 1. (40) and a f o r t i o r i

Pr{3ic1%1+:

X iE B ~ ~

x ~ = 1. (41)

This proves (ii) for the case that I and II are satisfied, so that it remains to establish the latter facts.

We recall that to prove l, we have to show that there is an e ( B ~ such that f o r all e<~ ~(B ~

inf P r { 3 i ~ : X ' E B ~ 1 7 6 (42)

x c ,at

C h o o s e any e > 0, and y, z z A~ with Ily - z II ~ e. Let

H,

be a h y p e r c u b e with centre at the origin and edges o f length

e / , / d

oriented along the c o o r d i n a t e axes, and let /-7 c H, be an arbitrary h y p e r c u b e whose edges are parallel to the c o o r d i n a t e axes as well. N o w construct the hyperrectangle

R,.:

oriented along the c o o r d i n a t e axes, which contains y and z as vertices. Then, since S is convex, and each point in A~

is at distance at least e from the b o u n d a r y o f S, it is easily verified that s + H~ = S for each vertex s (including y and z) o f

R,.:. The

edges o f

R,,-

are perpendicular.

Thus, if the vertices s ~ and

s 2

only differ in the j - t h coordinate, we can reach each point in

s2+H,,

from each point in

sJ+H,,

by a m o v e along the direction ei (j = 1 , . . . ,d). Hence, since the c o o r d i n a t e directions algorithm searches along one o f the c o o r d i n a t e vectors with equal probability, and since a next interior point is chosen uniformly on the linesegment c o n n e c t i n g the two previous hitpoints, we have that, for all x 6 y + H , ,

(12)

H.C.P. Berbee et aL / Hit-and-run algorithms to identify nonredundant constraints 195

/z(/q) (43)

Pr{X J e z + / q J X ~ = x}/> (drs)a

w h e r e rs is the m a x i m a l d i a m e t e r o f S (cf. d e f i n i t i o n (19)).

Next, let y a n d z be a r b i t r a r y p o i n t s in A~, not n e c e s s a r i l y s a t i s f y i n g NY- z[[ <~ e.

It is c l e a r that we can c h o o s e q & r , / e p o i n t s x ~ , . . . , x ~ on the straight line l c o n n e c t i n g y a n d z, such t h a t each two s u c c e s s i v e p o i n t s are at d i s t a n c e at m o s t e.

F u r t h e r m o r e , since S is c o n v e x , A , is c o n v e x , a n d /c= AF. Thus, a p p l y i n g (43) q + 1 times, we find that for all x c y + H , , a n d , for all Borel sets E such that p . ( E c v z + H ~ ) > O ,

Pr{X Iq+l~J c E IX ~ x}

> ~ P r { X ( q + ~ ) a c E c ~ z + H ~ , X ~ ' l c x i + H ~ ; i = ! , . . . . q ] X ~

with

>~6#(E ~ z + H~.)>O (44)

3 = ~ \ ~ ) (45)

Since the a b o v e l o w e r b o u n d d o e s not d e p e n d on x, it follows t h a t for all Borel sets E with ~ ( E ~ A ~ + H F ) > 0 that

i n f P r { 3 i ~ N * : X i ~ E I X ~ (46)

We now r e t u r n to o u r fixed Borel set B ~ with /x(B ~ > 0. It is i m m e d i a t e that there is an e ( B ~ such that, for all e <~ e(B~ p.(B'~ca A , + H,.) > 0. H e n c e , (46) is satisfied for E = B ~ w h i c h p r o v e s I.

F i n a l l y , we will p r o v e II, i.e. if 7/ is b o u n d e d f r o m a b o v e b y / x , then for o u r fixed a ~ t h e r e exists an e ( a ~ s u c h t h a t for all e <~ e(c~ ~

i ~ 0

Pr~{{X }i=oc A~ infinitely o f t e n } > 1 - o ~ . (47) Since p . ( S ) < o o , we c a n c h o o s e e ( d ~) such that Clx(S-Ao~#,~)<~t ~ Since tx is i n v a r i a n t this c h o i c e i m p l i e s that, if e ~< e(c~~

P r , { X ~ c S - A~} <~ c p . ( S - A, ) < c~ ~ (48) for all n. Thus, if

_N & sup{n IX" c A , } (49)

n>~d

t h e n , b y (48),

P r , 7 { N < n}<~ P r v { X " e S - A , } ~ a ~ (50) for all n >/d, so that

Pr,~{N = ~ } > 1 - o~ ~ (51)

w h i c h p r o v e s II.

(13)

196 H.(.'.P. Berbee et al. ! Hit-and-run algorilt m s to identi]~v n o n r e d u n d a n t constraints

4. A Bayesian stopping criterion

Let the n u m b e r of n o n r e d u n d a n t constraints of a system o f m linear inequalities be k~< m. We already observed that unless k = m and m distinct n o n r e d u n d a n t constraints have been f o u n d , it remains uncertain if all n o n r e d u n d a n t constraints have been identified. Let ~ be the n u m b e r o f hitpoints, and define

w

as the n u m b e r o f

distinct

n o n r e d u n d a n t constraints which have been f o u n d in the course o f these trials. (Notice that after n iterations o f an hit and run algorithm f i = 2 n ) . In this section we will develop a

Bayesian stopping criterion

which determines for each pair (& w) if the search p r o c e d u r e should be terminated, or not.

Assume that the k n o n r e d u n d a n t constraints are labeled 1 , . . . , k, and for each x c S ~ define (~(x) as the subset o f the hypersphere with centre x, such that a search from x in the direction o f _s will yield a hitpoint on the n o n r e d u n d a n t constraint with index i (i = 1 . . . k). D e n o t e the ( d - l ) - d i m e n s i o n a l Lebesgue measure by

L,l

~{ 9 }, and assume without loss o f generality that the ( d - l ) - d i m e n s i o n a l Lebes- gue measure of the surface area o f the above hypersphere is equal to 1. Then the probability that the ~-th search o f the

hypersphere directions algorithm

will yield a hitpoint on the i-th n o n r e d u n d a n t constraint is given by

0 ~ ( f i ) = f La

, { s 1 7 6 ~

( i = 1 . . . k). (52)

, ) 5

Next, define G , , c S as the set o f points in S from which the i-th n o n r e d u n d a n t constraint can be f o u n d by a search along direction es (i = 1 . . . k; j = 1 . . . . , d).

Then, a n a l o g o u s l y to the above reasoning, we obtain for the

coordinate directions algorithm

that the hitting probabilities are given by

O~(~)=~Tdl ~" Pr{X,,cG,,lXO=xO }

( i = 1 , . . . . , k ) . (53)

i - I

Hence, a p p l y i n g T h e o r e m s 2 and 3 to respectively (52) and (53) we obtain that given X" = x ~

pm 0,H(fi)= f L,,_~{~,(x)}p.(dx) ( i = 1 . . . . , k) (54)

and

1 a

lira 0 ~ ( ~ ) = - - x~ g ( G u ) ( i = 1 , k). (55)

,~-, 2d i~l . . . .

Hence, the hitting probabilities are asymptotically.lbced; the c o n v e r g e n c e rates are addressed in Smith (1984). For the two versions of the hit-and-run algorithm we denote these limiting values 0~ and

0~

respectively; observe that 0~ is not necessarily equal to

0~, (i = 1 . . . . , k).

In order to be able to p r o c e e d we will exploit this result by assuming that the actual hitting probability o f n o n r e d u n d a n t constraint i at each trial o f the h y p e r s p h e r e and c o o r d i n a t e directions algorithm is equal to the

(14)

H.C.P. Berbee et al./ Hit-and-run algorithms to identify nonredundant constraints 197 c o r r e s p o n d i n g a s y m p t o t i c probability 01 ~ and 0~ c respectively (i = 1 , . . . , k). Then, if the values of these probabilities would be given, several p r o p e r s t o p p i n g rules are obvious. In that case, a p r o c e d u r e could be terminated, for example, if the total probability o f the observed n o n r e d u n d a n t constraints exceeds a prescribed value, or if the total probability o f the observed n o n r e d u n d a n t constraints is equal to 1, i.e. if all k n o n r e d u n d a n t constraints have been found. The true n u m b e r o f n o n r e d u n - dant constraints of a system o f linear inequalities, and a f o r t i o r i the c o r r e s p o n d i n g hitting probabilities, are o f course frequently u n k n o w n . Therefore we a d o p t a statistical a p p r o a c h in which the data p r o d u c e d by an hit-and-run algorithm are used to gain i n f o r m a t i o n about their values. O u r starting point is the a s s u m p t i o n that the actual hitting probabilities at each trial are equal to the a s y m p t o t i c prob- abilities. Then the o u t p u t o f an algorithm is a sample from a multinomial distribution:

each cell o f the distribution c o r r e s p o n d s to a n o n r e d u n d a n t constraint, and the cell probabilities are equal to the c o r r e s p o n d i n g hitting probabilities. Thus, if we make no further notational distinction between 0, H and O~(i = 1 , . . . , k), the joint probabil- ity that in fi searches the ith n o n r e d u n d a n t constraint will be f o u n d n~ times ( i = 1 . . . . , k) is equal to

ti! k

p ( n , . . . . , n k l k , O, . . . 0k) l-[k |[ 0,, n (56)

z.,~=~0~ 1, v)=~n~ f i ) . G i v e n a s a m p l e ( n ~ , . . . , n ~ ) t h e m u l t i n o m i a l f o r m u l a ( 5 6 ) would enable us to learn a b o u t the values o f the u n k n o w n s (k, 0~, . . . , Ok). However, since it is u n k n o w n in a d v a n c e which o f the m constraints are n o n r e d u n d a n t , it is impossible to distinguish between different samples ( t h , . . . , nk) up to a relabeling o f the n o n r e d u n d a n t constraints. For example, if in fi = 5 searches one n o n r e d u n d a n t constraint has been f o u n d 4 times and a n o t h e r one once, it is u n k n o w n whether we observed ( n ~ , n , ) = ( 4 , 1 ) , (n~,n_~)=(1,4) or ( n l , n : , n 3 , n4, ns, n6, n v ) = (0, 4, 0, 1,0, 0, 0) etc. Thus we have to restrict ourselves to distinguishable aggregates o f the sample o u t c o m e s ( n ~ , . . . , nk) that are i n d e p e n d e n t o f the labeling o f the n o n r e d u n d a n t constraints and which do not contain n / s which are equal to 0. Given the result o f ~ searches we d e n o t e the a p p r o p r i a t e aggregates by { n ~ , . . . , n,,.} (recall that w is the n u m b e r o f observed n o n r e d u n d a n t constraints). The required probability o f these aggregates ~ n ~ , . . . , n,,} is given by the generalized multinomial distribution.

Theorem 4 ( B o e n d e r and R i n n o o y Kan, 1983a; Boender, 1984). Let a system ~1"

linear inequalities be given with k nonredundant constraints with probabilities 0~, . . . , Ok. Then the probability that in ii trials w d(Oerent nonredundant constraints are found, ~1 which one constraint is.lound n~ times, another constraint n, times etc., is given by the generalized multinomial distribution

1 f i ! "

- [1" E [1 o;:

p I { n , , . . . , nw)lk, o, . . . ok) I 1 7 : , ci! , = , n,!,• ... ,..,~.~-.1 , , (57)

(15)

198 H.C.P. Berbee et aL / Hit-and-run algorithms to ident(l~v nonredundant constraints

(V~ 10i = 1, Z~'~ln~ = r], n i > 0 ; i = 1 , . . . , w), where

cj ~= the number qf n,'s that are equal t o j (j = 1 . . . fi), (58) Sk [ w] ~- the set of all permutations of w d!fferent elements of

the set {l . . . . , k}. (59)

Now, (57) can be used in a Bayesian a p p r o a c h in which the u n k n o w n s k, 0~,.. , 0k are assumed to be themselves r a n d o m variables K, 0~ . . . . , 0 K for which a prior distribution can be specified. Given the result o f an hit-and-run algorithm, Bayes' rule is used to c o m p u t e the posterior distribution o f K, which i n c o r p o r a t e s both the prior beliefs and the sample information about the true n u m b e r o f n o n r e d u n d a n t constraints.

Theorem 5 ( B o e n d e r and R i n n o o y Kan, 1983b; Boender, 1984). Under the assumption of an arbitrary prior distribution p( 9 )./br the number of nonredundant constraints K, and, conditional on K = k, a symmetric Dirichlet prior with hyperparameter o~ for the hitting probabilities 0 ~ , . . . , Ok, i.e.

k

_ / ~ k - 1 ) ~ I1 01'-' (60)

p ( O l , . . . , O k ] k ) ( ( a _ 1)!)k i_i

the marginal posterior distribution of the true number of nonredundant constraints, conditional on an observed aggregate { n ~ , . . . , nw}, is equal to

p ( k [ { n , , . . . , n~.})oc p ( k ) ( a k - l ) ! k !

( f i + c ~ k - 1 ) ! ( k - w)! ' (61) where oc denotes proportionality.

O u r s t o p p i n g rule will be based on the k n o w l e d g e a b o u t the n u m b e r o f n o n r e d u n - dant constraints c o n t a i n e d in the posterior distribution (61). Observe that (61) is i n d e p e n d e n t o f the n u m b e r of times n~ that each o f the observed n o n r e d u n d a n t constraints has been f o u n d ( i = 1 , . . . , w), but only involves the total n u m b e r o f hitpoints ~ and the n u m b e r o f distinct observed n o n r e d u n d a n t constraints w, so that also o u r stopping rule only depends on the pair (ri, w).

Before we can apply (61) it remains to c h o o s e a p r o p e r prior distribution, Since a system o f m linear inequalities in a d dimensional space consists o f at least d + 1, and at most o f m n o n r e d u n d a n t constraints, the range o f the prior for K is d + 1 , . . . , m. The prior probability that the true n u m b e r o f n o n r e d u n d a n t constraints g is equal to k is a s s u m e d to grow linearly with k. Thus

p ( k ) o c k ( k = d + l . . . . , m ) . (62)

N o t e that in a Bayesian context a prior distribution is frequently c h o s e n uniform.

Since o u r analysis is based on the a s s u m p t i o n that the actual hitting probabilities

(16)

H. C.P. Berbee et aL / Hit-and-run algorithms to identify nonredundant constraints 1 9 9

are equal to the a s y m p t o t i c hitting probabilities, we deliberately chose a prior that will result in longer r u n n i n g times o f the algorithms than would be obtained for a u n i f o r m prior for K. H o w e v e r , the reader need not follow o u r suggestion and m a y c h o o s e any other version for (62), based on his own system o f linear inequalities.

Given K = k, we are only free to c h o o s e the h y p e r p a r a m e t e r c~ o f the symmetric Dirichlet prior distribution for the hitting probabilities 0 ~ , . . . , 0k (cf. (60)). The value o f a is a measure o f the deviation o f the prior (60) from its expectation

E ( O i ) = l / k

( i = l , . . ., k) (cf. Wi[ks, 1962). F o r a = l the s y m m e t r i c Dirichlet c o r r e s p o n d s to the u n i f o r m distribution on the unit simplex (~k=j 0, = 1; 0 < 0~ < 1;

i = 1 , . . . , k). For a > 1 a u n i q u e m a x i m u m is attained at the expected value, whose size increases if a gets larger: if a = oo then all probabilities are a priori assumed to be equal to

l / k

with probability 1: For 0 < a < 1 the distribution attains a u n i q u e m i n i m u m at the expected value which decreases as a a p p r o a c h e s 0. Preliminary diagnostic experiments s h o w e d that the posterior (61), and a f o r t i o r i o u r stopping rule, are sensitive to the choice for a. Hence, a user who is uncertain about the correct value for c~ is in a difficult position. To cope with this problem we run a hit-and-run algorithm a certain n u m b e r o f iterations and estimate a by G o o d ' s formula ( G o o d , 1965):

w )

i = l \ f f /

- l + w

i~l \ ff /

Then this estimate is used as if it is the true value c o r r e s p o n d i n g to the system o f linear inequalities u n d e r investigation.

N o w , given a choice o f the prior distribution and an observed sample o f hitpoints we can c o m p u t e the posterior distribution o f the true n u m b e r o f n o n r e d u n d a n t constraints f r o m (61). Then we can easily calculate the

posterior expected value

o f the n u m b e r o f n o n r e d u n d a n t constraints, which is well k n o w n to be the

optimal Bayesian estimate

with respect to a

quadratic loss function

(Lindley, 1978). For most (if, w) pairs, however, this m a y yield a real valued estimate, whereas the true n u m b e r o f n o n r e d u n d a n t constraints is evidently an integer. Therefore, since it easily shown that the optimal integer Bayesian estimate u n d e r a quadratic loss function is the round-off o f the real valued estimate, we will terminate our algorithms when this r o u n d - o f f o f the real valued optimal Bayesian estimate o f the n u m b e r o f n o n r e d u n - dant constraints is equal to the n u m b e r o f distinct n o n r e d u n d a n t constraints observed. That is, the algorithms are s t o p p e d when the current (fi, w) pair satisfies

kZ ( ~ , k - l ) ! k !

. . . .

I,l~,.wl

( f i + ~ k - l ) ! ( k - w ) ! < w + ~ "

E ( K [ ( r L w)) - k

" ( d k - 1 ) ! k T

Z..,

k ... {,1~,,,.I ( f i + 3 z k - l ) ! ( k - w ) !

( 6 4 )

(17)

2 0 0 H.C.P. Berbee et al. / Hit-and-run algorithms to identil]v nonredundant constraints

5. Computational results

In this section we describe the performance of the two hit-and-run algorithms and the Bayesian stopping rule on 3 practical and 3 randomly generated test problems. Note that the algorithms require an initial interior point. 11" not readily available, this may be obtained from a phase 1 procedure. This starting problem is"

a common burden of all nonredundancy identification procedures (cf. Karwan et al., 1983).

The most important features of our test problems are displayed in Table 1. The experimental design for the randomly generated test problems is taken from Karwan et al. (1983). The practical problems A and B are from Tischer (1968), practical problem C is from Meyerman (1966). The practical problems incorporate x~ 1> 0 ( i = 1 , . . . , d) as part of the set of constraints. The density of the problems (i.e., the fraction of non-zero coefficients of the constraints) is denoted by 6.

In Table 2 the number of seconds is shown that the hit-and-run algorithms required on a D E C 2060 computer to generate ~ = 100000 hitpoints; the table shows the estimates of the hyperparameter ~t as well. The Figures 3 up to 8 depict the evolution of the fraction of observed nonredundant constraints in the course of sampling:

stopping times (in seconds) are shown in Table 3. HD denotes the hypersphere directions algorithm and CD is the coordinate directions algorithm.

Figures 3-8 and Tables 1-3 reveal substantial differences between the hypersphere and coordinate directions algorithm. A final issue of theoretical and practical relevance is to what extent these algorithms are successful in approximating the uniform distribution over S in a limited number of experiments. Theoretically, exponential speed of convergence was established for the hypersphere directions method in Smith (1984), and we conjecture that a similar result is true for the coordinate directions method. Experimentally, we compared their performance on two 2-dimensional polytopes, of which the first one (Figure 9) was chosen to be very disadvantageous for the coordinate directions method and the other one (Figure 10) to be very advantageous. The restllts for 100 and 500 iterations of both methods are depicted in the figure. In spite of its simplicity, the coordinate directions method presents a very acceptable short run picture, providing additional practical confirma- tion of its superiority over its competitor.

6. Concluding remarks

Our experiments show that the hypersphere directions algorithm is inferior to its rival with respect to the required computer time to generate a given number of hitpoints as well as with respect to the number of identified nonredundant constraints per hitpoint. We will therefore in this section refer only to the coordinate directions algorithm. Although the best deterministic method does outperform this method on problems of small size (Karwan et al., 1983), it can be applied to (very) large

(18)

H.C.P. Berbee et aL / Hit-and-run algorithm.~ to ident!{i, nonredundant constraints 201

T a b l e 1 Test p r o b l e m s

P r o b l e m m d k ,,5

A 91 26 52 0.33

B 29 5 11 0.78

C 59 23 54 0.15

D 20 I0 19 1.00

E 20 10 18 0.50

F 30 10 18 0.50

T a b l e 2

C o m p u t a t i o n t i m e s f o r 1 0 0 0 0 0 h i t p o i n t s , a n d e s t i m a t e s f o r a

P r o b l e m A l g o r i t h m C o m p u t a t i o n d, time

A H D 1747 0.62

C D 61 4.18

B H D 315 1.55

C D 49 8.65

C H D 721 0.49

C D 20 1.02

D H D 577 0.71

C D 24 0.71

E H D 572 1.02

C D 24 0.79

F H D 813 1.01

C D 45 0.79

T a b l e 3 S t o p p i n g t i m e s

P r o b l e m A l g o r i t h m T i m e ~ w

A H D 1041 59 600 51

k = 52 C D 0.30 500 48

B H D 16.7 5300 7

k - I1 C D 0.05 100 11

C H D 721 100 000 49

k - 54 C D 1.1 5400 54

D H D 5.2 9 0 0 11

k - 19 C D 0.45 2 0 0 0 18

E H D 1.7 300 11

k - 18 C D 0.26 1100 15

F H D 2.4 300 11

k = 18 C D 0.5 1100 15

(19)

202 H . ( : P, [,lerbee et al./ Hit-and-run algorithm.s" to idemif[v nom'edundanl conslraints w / k

Hypersphere - -

t , Coordinate . . .

1.0

r).

i'

r

n.7

0 . 6

0 . 5 1000

/

i ] I I l l

3000 5000 70100 90100 11000 13000 Fig. 3. Problem A.

w l k

0.9[

0.

0, 0.

O. I _ _ ~ ! l I ,I I } l I

500 1500 2500 3500 4500 5500

rl

' 15ooo "~

n

6500 7 5 0 0 "

w l k ,I 1.0

0 . 9 -"

0 . 8

0 . 7

0 . 6

0 . 5

Fig. 4. Problem B.

n

4000 8000 12000 ]6000 20000 24000 28000 32000 36000 40000

/

/

i ~ 1 i i i

Fig. 5. Problem C.

(20)

H.C.P. Berbee et a l . / Hit-and-run algorithm.y to ident~/~r nonredundant constraints

'.,, / k H y p e r s p h e r e

~1 C o o r d i n a t e . . . .

1 , 0

0 , 9

0 , 8

0.7

0.6

I

0 . 5

j t j /

J J J J

/

J

I I L I I r I I I

/'000 8000 12000 16000 20000 24000 28000 32000 36000 40000

w/k

4

l.O

0.9

0.8

0.7

0.6

0.5

w/k A 1.0

0.9

0.8

0.7 /

0 . 6

0 . 5

Fig. 6. P r o b l e m D.

r . . .

/

I I I I t I f I 1

2500 5000 7500 10000 1250(t 15000 17500 20000 22500 25000

Fig. 7. P r o b l e m E.

n

v

n

n

I i t I I i i i f

v

4000 8000 12000 16000 20000 24000 28000 32000 36000 400UO

Fig. 8. P r o b l e m F.

203

(21)

204 H.C. P, Berhee et al./ Hit-and-run algorithms to identify nonredundant constraints

24

2 2

20

1 8

16

14

12

10

8 6

4

2

0

24

22

20

1 8

16

14

12

10

8

6

4 2

0

u . . . s p h e r e d i r e c t i o n s ,ints

24 ~ - o r d i n a t e d i r e c t i o n s 3 p o i n t s

22

18

16

14

12

10 20

6

4

2

I I I I I I I 0

0 2 4 6 8 10 12 14 16 18

H y p e r s p e r e directions

~ 500 points

24

22

20

18

16

1 ~,

12

10

8

,.:: .~

~'~ 2

I I I I ~ l I I 0

2 4 6 8 10 12 14 16 18

i I - I I I I r I

0 2 4 6 8 10 12 14 16 18

C o o r d i n a t e directions 500 points

;'L'J'~

~ , ,..

~".e

;* % . . , -

::.::

:~yK.

P "'1

i i I i i I I

4 6 8 10 12 14 16 18

Fig. 9.

0021

(22)

H.C.P. Berbee et al. / Hit-and-run algorithms to ident(~v nonredundant constraints 205

24

22

2 0

18

1 6

1 4 -

1 2 -

1 0 -

8 -

6 -

4 -

2

0 ,

24

22

20

18

16

14

12

10

8

6

4

2

0

H y p e r s p h e r e d i r e c t i o n s 100 p o i n t s

89 4 ~ ~ Ib 12 14 16 18 0

H y p e r s p h e r e d i r e c t i o n s 500 p o i n t s

24

22

20

18

16

14

12

10

8

6

4

24- C o o r d i n a t e d i r e c t i o n s 100 p o i n t s

22-

20-

18-

16-

l g -

10

8 6

2

2 4 6 1 12 14 16 18

C o o r d i n a t e directions 500 points

2 ;

i i I 1 I ] I

0 2 4 6 8 10 12 14 16 18

I i I I I i i I I 0

2 4 6 8 10 12 14 16 18

Fig. I0.

(23)

206 H.C.P. Berbee et al. / Hit-and-run algorithms to ident!fy nonredundant constraints

problems and therefore was judged in Karwan et al. (11983) to be one o f the most attractive practical possibilities to test for redundancy.

So far the random method has been viewed as a preprocessor that is capable o f identifying nonredundant constraints. In an optimization algorithm for (possibly nonlinear) objective functions with linear inequality constraints, the nonidentified constraints can be omitted. Obviously, the solution produced by the optimization algorithm has to be checked for feasibility afterwards, since the omission o f s o m e constraints may not be justified. In case o f linear programming, feasibility can be easily restored by dual simplex steps.

It is well k n o w n that the optimal solution o f a linear programming problem is attained in one o f the vertices o f the feasible region. Therefore for linear programming applications one may investigate the possibility of modifying the coordinate direc- tions algorithm to m o v e in the direction o f the optimal vertex. The modified method may even converge rapidly to this optimal vertex. If not, then at least a better preprocessor will be obtained, since the nonredundant constraints which are active in the optimal vertex will have relatively greater hitting probabilities.

Finally, for nonlinear objective functions with linear inequality constraints the optimal solution may be situated anywhere in the feasible region. Since the coordin- ate directions algorithm involves an extremely fast method to generate (asymptoti- cally) uniform points over such a feasible region, the method offers a starting point for various constrained

global optimization

procedures in which the generation o f such points w o u l d be a first step (Timmer, 1984).

References

R.B. Ash, Real Ana(vsis and Probability (Academic Press, New York, 1972).

M.L. Balinski, "'An algorithm for finding all vertices of convex polyhedral sets," Journal of the Socie O, q]" Industrial Applied Mathematics 9 ( 1961 ) 72-88.

C.G.E. Boender and A.H.G. Rinnooy Kan, "'A Bayesian analysis of the number of cells of a multinomial distribution,'" The Statistician 32 (1983a) 24(I-248.

C.G.E. Boender and A.H.G. Rinnooy Kan, "Bayesian estimation of animal population size," Report 8322/0 Econometric Institute, Erasmus University Rotterdam (Rotterdam, 1983b).

C.G.E. Boender, "'The generalized multinomial distribution: A Bayesian analysis and applications,"

Ph.D. Thesis, Erasmus University Rotterdam (Rotterdam, 1984).

A. Boneh and A. Golan, "'Constraints redundancy and feasible region boundedness by random feasible points generated," paper presented at EURO 111 (1979).

A. Boneh, "'A probabilistic algorithm identifying redundancy by a random feasible point generator (RFPG)," in: Karwan et al. (1983).

G.H. Bradley, G.G. Brown and G.W. Graves, "Structural redundancy in large scale optimization models,"

Technical Report CA 93940, Naval Postgraduate School ( Monterey, 1980); also in: Karwan et al. (1983).

A.L. Brearly, G. Mitra and H.P. Williams, "Analysis of mathematical programming problems prior to applying the simplex algorithm," Mathematical Programming 8 (1975) 54-83.

T. Gal, "Zur ldentifikation redundanter Nebenbedingungen in linearen Programmen,'" Zeitschrift fiir Operations Research 19 (1975) 19-28.

I.J. Good, The Estimation o1" Probabilities, Research Monograph no. 30, (The M.I.T. Press, Cambridge, MA, 1965).

M.H. Karwan, V. Lotfi, J. Telgen and S. Zionts, eds., Redundancy in Mathematical Programming (Springer-Verlag, Berlin, 1983).

Figure

Updating...

References

Related subjects :