Dear All,
Thank you for the replies to my first thread here:
http://r.789695.n4.nabble.com/global-optimisation-with-inequality-constraints-td3799258.html.
So far the best result is achieved via a penalised objective function. This
was suggested by someone on this list privately. I am still looking into some
of the options mentioned in the original thread, but I have been advised that
there may be better ways if I present the actual problem with a reproducible
example.
In principle the problem can be solved by linear programming, so I include code
for my attempt at this via RGLPK. It says that there is no feasible solution,
but the solution is known analytically in the case below.
Here is the precise problem:
Minimise, over 100×1 real vectors v,
Max_i(|v_i|) such that X'v=e_2,
where X is a given 100×2 matrix and e_2 =(0,1)'. The v_i are the elements of v.
I have put the actual X matrix at the end of this post, along with a feasible
starting value for v.
The correct minimum is 0.01287957, obtained with v_i=0.01287957 for i<=50 and
v_i = 0.01287957 for i>=51.
Here is the code for the penalised objective function approach, adapted very
slightly from what someone on this list sent me:
......................................................................
X <- # See end of this message for the X data
x1 <- X[, 1]
x2 <- X[, 2]
fun <- function(q) {
mu <- 0.1
max(abs(v)) + (sum(v*x1)^2 + (1-sum(x2*v))^2)/(2*mu)
}
vstart <- # feasible starting value. See end of this post.
sol <- optim(vstart, fun, method="L-BFGS-B", lower=rep(-1, 100),
upper=rep(1,100))
max(abs(sol$par))
.........................................................................
This gets quite near, around 0.015-0.016 for me by varying mu.
Alternatively the problem can be set up as a linear programming task as follows:
Minimise, over 100×1 real vectors v, and over scalars q >= 0,
q
such that, for i=1,...,100
v_i<=q
v_i>=-q
X'v=e_2
Here is my RGLPK code:
.........................................................................
X <- # See end of this message for the X data
XROWS <- 100
XCOLS <- 2
e_2=rep(0,times=XCOLS)
e2[2]<- 1
obj <- c(rep(0,XROWS),1) # coefficients on v_1, . . . , v_100, q.
mat <- matrix(rbind(cbind(diag(XROWS), rep(-1,XROWS)), cbind(diag(XROWS),
rep(1,XROWS)), cbind(t(X), rep(0,XCOLS)), cbind(t(rep(0,XROWS)), 1)),
nrow=2*XROWS+XCOLS+1)
dir <- c(rep("<=", XROWS), rep(">=", XROWS), rep("==", XCOLS), ">=")
rhs <- c(rep(0, 2*XROWS), e_2, 0)
sol <- Rglpk_solve_LP(obj, mat, dir, rhs, types = NULL, max = FALSE,
bounds = c(-5,5), verbose = TRUE)
...........................................................................
The output is
"
GLPK Simplex Optimizer, v4.42
203 rows, 101 columns, 601 non-zeros
0: obj = 0.000000000e+000 infeas = 1.000e+000 (2)
4: obj = 0.000000000e+000 infeas = 1.000e+000 (1)
PROBLEM HAS NO FEASIBLE SOLUTION
"
I have also tried setting the problem up with a small interval around the
equality constraints rather than having strict equalities, but could not get
the correct solution this way either.
Maybe I am making an error with RGLPK - I have been told it should work for a
problem of this size and much larger. I have also tried DEoptim, IpSolve and
ConstrOptim.
Regards,
Gareth
...............................................................................
Values of X and vstart below
...............................................................................
vstart=
-0.025251183
-0.022301089
-0.020429759
-0.01902228
-0.017877586
-0.016903415
-0.016049376
-0.015284788
-0.014589517
-0.0139496
-0.01335494
-0.012797983
-0.012272922
-0.011775194
-0.011301139
-0.010847778
-0.010412649
-0.009993692
-0.009589163
-0.009197573
-0.00881764
-0.008448248
-0.008088421
-0.007737298
-0.007394116
-0.007058194
-0.006728922
-0.006405748
-0.006088173
-0.005775744
-0.005468043
-0.005164689
-0.00486533
-0.004569638
-0.00427731
-0.003988063
-0.003701629
-0.003417759
-0.003136215
-0.002856772
-0.002579217
-0.002303343
-0.002028954
-0.00175586
-0.001483876
-0.001212825
-0.00094253
-0.00067282
-0.000403526
-0.000134481
0.000134481
0.000403526
0.00067282
0.00094253
0.001212825
0.001483876
0.00175586
0.002028954
0.002303343
0.002579217
0.002856772
0.003136215
0.003417759
0.003701629
0.003988063
0.00427731
0.004569638
0.00486533
0.005164689
0.005468043
0.005775744
0.006088173
0.006405748
0.006728922
0.007058194
0.007394116
0.007737298
0.008088421
0.008448248
0.00881764
0.009197573
0.009589163
0.009993692
0.010412649
0.010847778
0.011301139
0.011775194
0.012272922
0.012797983
0.01335494
0.0139496
0.014589517
0.015284788
0.016049376
0.016903415
0.017877586
0.01902228
0.020429759
0.022301089
0.025251183
.....................................................................................................................................
X=
1 -2.330078923
1 -2.057855981
1 -1.885177032
1 -1.755300501
1 -1.649672679
1 -1.559779992
1 -1.480972651
1 -1.410419531
1 -1.346262665
1 -1.287213733
1 -1.232340861
1 -1.180947041
1 -1.13249653
1 -1.086568115
1 -1.042824239
1 -1.000989917
1 -0.960837931
1 -0.922178178
1 -0.884849841
1 -0.848715527
1 -0.813656808
1 -0.779570774
1 -0.746367337
1 -0.713967098
1 -0.682299633
1 -0.651302112
1 -0.62091817
1 -0.591096977
1 -0.561792466
1 -0.532962693
1 -0.504569287
1 -0.476576998
1 -0.448953298
1 -0.421668052
1 -0.39469322
1 -0.368002611
1 -0.341571661
1 -0.315377237
1 -0.289397474
1 -0.263611615
1 -0.237999879
1 -0.212543342
1 -0.187223821
1 -0.162023779
1 -0.136926226
1 -0.111914639
1 -0.08697288
1 -0.062085116
1 -0.037235755
1 -0.012409369
1 0.012409369
1 0.037235755
1 0.062085116
1 0.08697288
1 0.111914639
1 0.136926226
1 0.162023779
1 0.187223821
1 0.212543342
1 0.237999879
1 0.263611615
1 0.289397474
1 0.315377237
1 0.341571661
1 0.368002611
1 0.39469322
1 0.421668052
1 0.448953298
1 0.476576998
1 0.504569287
1 0.532962693
1 0.561792466
1 0.591096977
1 0.62091817
1 0.651302112
1 0.682299633
1 0.713967098
1 0.746367337
1 0.779570774
1 0.813656808
1 0.848715527
1 0.884849841
1 0.922178178
1 0.960837931
1 1.000989917
1 1.042824239
1 1.086568115
1 1.13249653
1 1.180947041
1 1.232340861
1 1.287213733
1 1.346262665
1 1.410419531
1 1.480972651
1 1.559779992
1 1.649672679
1 1.755300501
1 1.885177032
1 2.057855981
1 2.330078923
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
