1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
|
Rcgminb <- function(par, fn, gr, lower, upper, bdmsk = NULL, control = list(), ...) {
## An R version of the conjugate gradient minimization
## using the Dai-Yuan ideas
#
# Input:
# par = a vector containing the starting point
# fn = objective function (assumed to be sufficeintly
# differentiable)
# gr = gradient of objective function
# lower = vector of lower bounds on parameters
# upper = vector of upper bounds on parameters
# Note: free parameters outside bounds will be adjusted to
# bounds.
# bdmsk = control vector for bounds and masks. Parameters
# for which bdmsk are 1
# are unconstrained or 'free', those with bdmsk 0 are
# masked i.e., fixed.
# For historical reasons, we use the same array as an
# indicator that a
# parameter is at a lower bound (-3) or upper bound (-1)
# control = list of control parameters
# maxit = a limit on the number of iterations (default 500)
# maximize = TRUE to maximize the function (default FALSE)
# trace = 0 (default) for no output,
# >0 for output (bigger => more output)
# eps=1.0e-7 (default) for use in computing numerical
# gradient approximations.
# dowarn=TRUE by default. Set FALSE to suppress warnings.
#
# Output:
# A list with components:
#
# par: The best set of parameters found.
#
# value: The value of 'fn' corresponding to 'par'.
#
# counts: A two-element integer vector giving the number of
# calls to
# 'fn' and 'gr' respectively. This excludes those calls
# needed
# to compute the Hessian, if requested, and any calls to
# 'fn'
# to compute a finite-difference approximation to the
# gradient.
#
# convergence: An integer code. '0' indicates successful
# convergence.
# Error codes are
# '0' converged
# '1' indicates that the function evaluation count
# 'maxfeval'
# was reached.
# '2' indicates initial point is infeasible
#
# message: A character string giving any additional
# information returned
# by the optimizer, or 'NULL'.
#
# bdmsk: Returned index describing the status of bounds and
# masks at the
# proposed solution. Parameters for which bdmsk are 1 are
# unconstrained
# or 'free', those with bdmsk 0 are masked i.e., fixed. For
# historical
# reasons, we indicate a parameter is at a lower bound
# using -3
# or upper bound using -1.
#
#
# Author: John C Nash
# Date: April 2, 2009; revised July 28, 2009
#################################################################
# control defaults -- idea from spg
ctrl <- list(maxit = 500, maximize = FALSE, trace = 0, eps = 1e-07,
dowarn = TRUE, tol=0)
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
npar<-length(par)
if (ctrl$tol == 0) tol <- npar * (npar * .Machine$double.eps) # for gradient test. Note -- integer overflow if n*n*d.eps
else tol<-ctrl$tol
maxit <- ctrl$maxit # limit on function evaluations
maximize <- ctrl$maximize # TRUE to maximize the function
trace <- ctrl$trace # 0 for no output, >0 for output (bigger => more output)
if (trace > 2) cat("trace = ", trace, "\n")
eps <- ctrl$eps
fargs <- list(...) # the ... arguments that are extra function / gradient data
grNULL <- is.null(gr)
dowarn <- ctrl$dowarn #
#############################################
if (maximize) {
warning("Rcgmin no longer supports maximize 111121 -- see documentation")
msg<-"Rcgmin no longer supports maximize 111121"
ans <- list(par, NA, c(0, 0), 9999, msg, bdmsk)
return(ans)
}
#############################################
# gr MUST be provided
if (is.null(gr)) { # if gr function is not provided STOP (Rvmmin has definition)
stop("A gradient calculation (analytic or numerical) MUST be provided for Rcgmin")
}
if ( is.character(gr) ) {
# Convert string to function call, assuming it is a numerical gradient function
mygr<-function(par=par, userfn=fn, ...){
do.call(gr, list(par, userfn, ...))
}
} else { mygr<-gr }
############# end test gr ####################
## Set working parameters (See CNM Alg 22)
if (trace > 0) {
cat("Rcgmin -- J C Nash 2009 - bounds constraint version of new CG\n")
cat("an R implementation of Alg 22 with Yuan/Dai modification\n")
}
bvec <- par # copy the parameter vector
n <- length(bvec) # number of elements in par vector
maxfeval <- round(sqrt(n + 1) * maxit) # change 091219
ig <- 0 # count gradient evaluations
ifn <- 1 # count function evaluations (we always make 1 try below)
stepredn <- 0.15 # Step reduction in line search
acctol <- 1e-04 # acceptable point tolerance
reltest <- 100 # relative equality test
ceps <- .Machine$double.eps * reltest
accpoint <- as.logical(FALSE) # so far do not have an acceptable point
cyclimit <- min(2.5 * n, 10 + sqrt(n)) #!! upper bound on when we restart CG cycle
fargs <- list(...) # function arguments
if (trace > 2) {
cat("Extra function arguments:")
print(fargs)
}
# set default masks if not defined
if (is.null(bdmsk)) {
bdmsk <- rep(1, n)
}
if (trace > 2) {
cat("bdmsk:")
print(bdmsk)
}
# Routine should NOT be called directly without bounds.
# Still do checks to get nolower, noupper, bounds
if (is.null(lower) || !any(is.finite(lower)))
nolower = TRUE
else nolower = FALSE
if (is.null(upper) || !any(is.finite(upper)))
noupper = TRUE
else noupper = FALSE
if (nolower && noupper && all(bdmsk == 1))
bounds = FALSE
else bounds = TRUE
if (trace > 2)
cat("Bounds: nolower = ", nolower, " noupper = ", noupper,
" bounds = ", bounds, "\n")
if (nolower)
lower <- rep(-Inf, n)
if (noupper)
upper <- rep(Inf, n)
######## check bounds and masks #############
## NOTE: do this inline to avoid call to external routine
if (bounds) {
# Make sure to expand lower and upper
if (!nolower & (length(lower) < n))
{
## tmp<-readline('Check length lower ')
if (length(lower) == 1) {
lower <- rep(lower, n)
}
else {
stop("1<length(lower)<n")
}
} # else lower OK
if (!noupper & (length(upper) < n))
{
if (length(upper) == 1) {
upper <- rep(upper, n)
}
else {
stop("1<length(upper)<n")
}
} # else upper OK
# At this point, we have full bounds in play
# This implementation as a loop, but try later to vectorize
for (i in 1:n) {
# cat('i = ',i,'\n')
if (bdmsk[i] == 0) {
# NOTE: we do not change masked parameters, even if out of bounds
if (!nolower) {
if (bvec[i] < lower[i]) {
wmsg <- paste(bvec[i], " = MASKED x[", i,
"] < lower bound = ", lower[i], sep = "")
if (dowarn)
warning(wmsg)
}
}
if (!noupper) {
if (bvec[i] > upper[i]) {
wmsg <- paste(bvec[i], " = MASKED x[", i,
"] > upper bound = ", upper[i], sep = "")
if (dowarn)
warning(wmsg)
}
}
}
else {
# not masked, so must be free or active constraint
if (!nolower) {
if (bvec[i] <= lower[i]) {
# changed 090814 to ensure bdmsk is set
wmsg <- paste("x[", i, "], set ", bvec[i],
" to lower bound = ", lower[i], sep = "")
if (dowarn && (bvec[i] != lower[i]))
warning(wmsg)
bvec[i] <- lower[i]
bdmsk[i] <- -3 # active lower bound
}
}
if (!noupper) {
if (bvec[i] >= upper[i]) {
# changed 090814 to ensure bdmsk is set
wmsg <- paste("x[", i, "], set ", bvec[i],
" to upper bound = ", upper[i], sep = "")
if (dowarn && (bvec[i] != upper[i]))
warning(wmsg)
bvec[i] <- upper[i]
bdmsk[i] <- -1 # active upper bound
}
}
} # end not masked
} # end loop for bound/mask check
} else stop("Do not call Rcgminb without bounds")
############## end bounds check #############
# Initial function value -- may NOT be at initial point
# specified by user.
if (trace > 2) {
cat("Try function at initial point:")
print(bvec)
}
f <- try(fn(bvec, ...), silent = TRUE) # Compute the function at initial point.
if (trace > 0) {
cat("Initial function value=", f, "\n")
}
if (inherits(f,"try-error")) {
msg <- "Initial point is infeasible."
if (trace > 0)
cat(msg, "\n")
ans <- list(par, NA, c(ifn, 0), 2, msg, bdmsk)
names(ans) <- c("par", "value", "counts", "convergence",
"message", "bdmsk")
return(ans)
}
fmin <- f
if (trace > 0)
cat("Initial fn=", f, "\n")
if (trace > 2)
print(bvec)
# Start the minimization process
keepgoing <- TRUE
msg <- "not finished" # in case we exit somehow
oldstep <- 0.8 #!! 2/3 #!!?? WHY?
####################################################################
fdiff <- NA # initially no decrease
cycle <- 0 # !! cycle loop counter
while (keepgoing) {
# main loop -- must remember to break out of it!!
t <- as.vector(rep(0, n)) # zero step vector
c <- t # zero 'last' gradient
while (keepgoing && (cycle < cyclimit)) {
## cycle loop
cycle <- cycle + 1
if (trace > 0)
cat(ifn, " ", ig, " ", cycle, " ", fmin, " last decrease=",
fdiff, "\n")
if (trace > 2) {
print(bvec)
cat("\n")
}
if (ifn > maxfeval) {
msg <- paste("Too many function evaluations (> ",
maxfeval, ") ", sep = "")
if (trace > 0)
cat(msg, "\n")
ans <- list(par, fmin, c(ifn, ig), 1, msg, bdmsk) # 1 indicates not converged in function limit
names(ans) <- c("par", "value", "counts", "convergence",
"message", "bdmsk")
return(ans)
}
par <- bvec # save best parameters
ig <- ig + 1
if (ig > maxit) {
msg <- paste("Too many gradient evaluations (> ",
maxit, ") ", sep = "")
if (trace > 0)
cat(msg, "\n")
ans <- list(par, fmin, c(ifn, ig), 1, msg, bdmsk) # 1 indicates not converged in function or gradient limit
names(ans) <- c("par", "value", "counts", "convergence",
"message", "bdmsk")
return(ans)
}
g <- mygr(bvec, ...)
if (bounds)
{
## Bounds and masks adjustment of gradient ##
## first try with looping -- later try to vectorize
if (trace > 2) {
cat("bdmsk:")
print(bdmsk)
}
for (i in 1:n) {
if ((bdmsk[i] == 0)) {
# masked, so gradient component is zero
g[i] <- 0
}
else {
if (bdmsk[i] == 1) {
if (trace > 1)
cat("Parameter ", i, " is free\n")
}
else {
if ((bdmsk[i] + 2) * g[i] < 0) {
# test for -ve gradient at upper bound, +ve at lower bound
g[i] <- 0 # in which case active mask or constraint and zero gradient component
}
else {
bdmsk[i] <- 1 # freeing parameter i
if (trace > 1)
cat("freeing parameter ", i, "\n")
}
}
}
} # end masking loop on i
if (trace > 2) {
cat("bdmsk adj:\n")
print(bdmsk)
cat("proj-g:\n")
print(g)
}
} # end if bounds
## end bounds and masks adjustment of gradient
g1 <- sum(g * (g - c)) # gradient * grad-difference
g2 <- sum(t * (g - c)) # oldsearch * grad-difference
gradsqr <- sum(g * g)
if (trace > 1) {
cat("Gradsqr = ", gradsqr, " g1, g2 ", g1, " ",
g2, " fmin=", fmin, "\n")
}
c <- g # save last gradient
g3 <- 1 # !! Default to 1 to ensure it is defined -- t==0 on first cycle
if (gradsqr > tol * (abs(fmin) + reltest)) {
if (g2 > 0) {
betaDY <- gradsqr/g2
betaHS <- g1/g2
g3 <- max(0, min(betaHS, betaDY)) # g3 is our new 'beta' !! Dai/Yuan 2001, (4.2)
}
}
else {
msg <- paste("Very small gradient -- gradsqr =",
gradsqr, sep = " ")
if (trace > 0)
cat(msg, "\n")
keepgoing <- FALSE # done loops -- should we break ??
break # to leave inner loop
}
if (trace > 2)
cat("Betak = g3 = ", g3, "\n")
if (g3 == 0 || cycle >= cyclimit) {
# we are resetting to gradient in this case
if (trace > 0) {
if (cycle < cyclimit)
cat("Yuan/Dai cycle reset\n")
else cat("Cycle limit reached -- reset\n")
}
fdiff <- NA
cycle <- 0
break #!!
}
else {
# drop through if not Yuan/Dai cycle reset
t <- t * g3 - g # t starts at zero, later is step vector
gradproj <- sum(t * g) # gradient projection
if (trace > 1)
cat("Gradproj =", gradproj, "\n")
if (bounds)
{
## Adjust search direction for masks
if (trace > 2) {
cat("t:\n")
print(t)
}
t[which(bdmsk <= 0)] <- 0 # apply mask constraint
if (trace > 2) {
cat("adj-t:\n")
print(t)
}
## end adjust search direction for masks
} # end if bounds
# ?? Why do we not check gradproj size??
########################################################
#### Line search ####
OKpoint <- FALSE
if (trace > 2)
cat("Start linesearch with oldstep=", oldstep,
"\n")
steplength <- oldstep * 1.5 #!! try a bit bigger
f <- fmin
changed <- TRUE # Need to set so loop will start
while ((f >= fmin) && changed) {
if (bounds)
{
# Box constraint -- adjust step length
for (i in 1:n) {
# loop on parameters -- vectorize??
if ((bdmsk[i] == 1) && (t[i] != 0))
{
# only concerned with free parameters and non-zero search
# dimension
if (t[i] < 0) {
# going down. Look at lower bound
trystep <- (lower[i] - par[i])/t[i] # t[i] < 0 so this is positive
}
else {
# going up, check upper bound
trystep <- (upper[i] - par[i])/t[i] # t[i] > 0 so this is positive
}
if (trace > 2)
cat("steplength, trystep:", steplength,
trystep, "\n")
steplength <- min(steplength, trystep) # reduce as necessary
} # end steplength reduction
} # end loop on i to reduce step length
if (trace > 1)
cat("reset steplegth=", steplength, "\n")
# end box constraint adjustment of step length
} # end if bounds
bvec <- par + steplength * t
changed <- (!identical((bvec + reltest), (par + reltest)))
if (changed) {
# compute newstep, if possible
f <- fn(bvec, ...) # Because we need the value for linesearch, don't use try()
# instead preferring to fail out, which will hopefully be
# unlikely.
ifn <- ifn + 1
if (is.na(f) || (!is.finite(f))) {
warning("Rcgmin - undefined function")
f <- .Machine$double.xmax
}
if (f < fmin) {
f1 <- f # Hold onto value
}
else {
savestep<-steplength
steplength <- steplength * stepredn
if (steplength >=savestep) changed<-FALSE
if (trace > 0)
cat("*")
}
}
} # end while
changed1 <- changed # Change in parameters occured in step reduction
if (changed1)
{
## ?? should we check for reduction? or is this done in if
# (newstep >0) ?
newstep <- 2 * (f - fmin - gradproj * steplength) # JN 081219 change
if (newstep > 0) {
newstep = -(gradproj * steplength * steplength/newstep)
}
if (bounds)
{
# Box constraint -- adjust step length
for (i in 1:n) {
# loop on parameters -- vectorize??
if ((bdmsk[i] == 1) && (t[i] != 0))
{
# only concerned with free parameters and non-zero search dimension
if (t[i] < 0) {
# going down. Look at lower bound
trystep <- (lower[i] - par[i])/t[i] # t[i] < 0 so this is positive
}
else {
# going up, check upper bound
trystep <- (upper[i] - par[i])/t[i] # t[i] > 0 so this is positive
}
if (trace > 2)
cat("newstep, trystep:", newstep,
trystep, "\n")
newstep <- min(newstep, trystep) # reduce as necessary
} # end newstep reduction
} # end loop on i to reduce step length
if (trace > 2)
cat("reset newstep=", newstep, "\n")
# end box constraint adjustment of step length
} # end if bounds
bvec <- par + newstep * t
changed <- (!identical((bvec + reltest),
(par + reltest)))
if (changed) {
f <- fn(bvec, ...)
ifn <- ifn + 1
}
if (trace > 2)
cat("fmin, f1, f: ", fmin, f1, f, "\n")
if (f < min(fmin, f1)) {
# success
OKpoint <- TRUE
accpoint <- (f <= fmin + gradproj * newstep *
acctol)
fdiff <- (fmin - f) # check decrease
fmin <- f
oldstep <- newstep # !! save it
}
else {
if (f1 < fmin) {
bvec <- par + steplength * t # reset best point
accpoint <- (f1 <= fmin + gradproj *
steplength * acctol)
OKpoint <- TRUE # Because f1 < fmin
fdiff <- (fmin - f1) # check decrease
fmin <- f1
oldstep <- steplength #!! save it
}
else {
# no reduction
fdiff <- NA
accpoint <- FALSE
} # f1<?fmin
} # f < min(f1, fmin)
if (trace > 1)
cat("accpoint = ", accpoint, " OKpoint = ",
OKpoint, "\n")
if (!accpoint) {
msg <- "No acceptable point -- exit loop"
if (trace > 0)
cat("\n", msg, "\n")
keepgoing <- FALSE
break #!!
}
} # changed1
else {
# not changed on step redn
if (cycle == 1) {
msg <- " Converged -- no progress on new CG cycle"
if (trace > 0)
cat("\n", msg, "\n")
keekpgoing <- FALSE
break #!!
}
} # end else
} # end of test on Yuan/Dai condition
#### End line search ####
if (bounds)
{
## Reactivate constraints?? -- should check for infinite
# bounds
for (i in 1:n) {
if (bdmsk[i] == 1)
{
# only interested in free parameters
if (is.finite(lower[i])) {
# JN091020 -- need to use abs in case bounds negative
if ((bvec[i] - lower[i]) < ceps * (abs(lower[i]) +
1))
{
# are we near or lower than lower bd
if (trace > 2)
cat("(re)activate lower bd ",
i, " at ", lower[i], "\n")
bdmsk[i] <- -3
} # end lower bd reactivate
}
if (is.finite(upper[i])) {
# JN091020 -- need to use abs in case bounds negative
if ((upper[i] - bvec[i]) < ceps * (abs(upper[i]) +
1))
{
# are we near or above upper bd
if (trace > 2)
cat("(re)activate upper bd ",
i, " at ", upper[i], "\n")
bdmsk[i] <- -1
} # end lower bd reactivate
}
} # end test on free params
} # end reactivate constraints
} # end if bounds
} # end of inner loop (cycle)
if (oldstep < acctol) {
oldstep <- acctol
}
# steplength
if (oldstep > 1) {
oldstep <- 1
}
if (trace > 1)
cat("End inner loop, cycle =", cycle, "\n")
} # end of outer loop
msg <- "Rcgmin seems to have converged"
if (trace > 0)
cat(msg, "\n")
# par: The best set of parameters found.
# value: The value of 'fn' corresponding to 'par'.
# counts: number of calls to 'fn' and 'gr' (2 elements)
# convergence: An integer code. '0' indicates successful
# convergence.
# message: A character string or 'NULL'.
# if (maximize)
# fmin <- -fmin
ans <- list(par, fmin, c(ifn, ig), 0, msg, bdmsk)
names(ans) <- c("par", "value", "counts", "convergence",
"message", "bdmsk")
return(ans)
} ## end of Rcgmin
|
