R Under development (unstable) (2024-10-16 r87241 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > if(!require("GNE"))stop("this test requires package GNE.") Loading required package: GNE Loading required package: alabama Loading required package: numDeriv Loading required package: nleqslv Loading required package: BB Loading required package: SQUAREM > > > > #------------------------------------------------------------------------------- > # (1) Example 5 of von Facchinei et al. (2007) > #------------------------------------------------------------------------------- > > dimx <- c(1, 1) > #O_i(x) > obj <- function(x, i) + { + if(i == 1) + res <- (x[1]-1)^2 + if(i == 2) + res <- 2*(x[2]-1/2)^2 + res + } > #Gr_x_j O_i(x) > grobj <- function(x, i, j) + { + if(i == 1) + res <- c(2*(x[1]-1), 0) + if(i == 2) + res <- c(0, 2*(x[2]-1/2)) + res[j] + } > #Gr_x_k Gr_x_j O_i(x) > heobj <- function(x, i, j, k) + 2 * (i == j && j == k) > > dimlam <- c(1, 1) > #constraint function g_i(x) > h <- function(x) + sum(x[1:2]) - 1 > > #Gr_x_j h(x) > jach <- function(x) + matrix(c(1, 1), 1) > > > x0 <- rexp(sum(dimx)) > y <- rexp(sum(dimx)) > xy1 <- c(y[1], x0[2]) > xy2 <- c(x0[1], y[2]) > alpha <- 0.1 > > resgap <- gapNIR(x0, y, dimx, obj=obj, echo=TRUE) > gapcheck <- obj(x0, 1) - obj(c(y[1], x0[2]), 1) - alpha/2*(x0[1]-y[1])^2 + obj(x0, 2) - obj(c(x0[1], y[2]), 2) - alpha/2*(x0[2]-y[2])^2 > > tol <- .Machine$double.eps^(1/2) > if(abs(resgap - gapcheck) > tol) + stop("wrong result 1") > > resgrgap <- gradxgapNIR(x0, y, dimx, grobj=grobj) > > grcheck <- c(grobj(x0, 1, 1) - grobj(xy1, 1, 1) + grobj(x0, 2, 1) - grobj(xy2, 2, 1), + grobj(x0, 1, 2) - grobj(xy1, 1, 2) + grobj(x0, 2, 2) - grobj(xy2, 2, 2)) > grcheck <- grcheck + c(grobj(xy1, 1, 1), grobj(xy2, 2, 2))- alpha*(x0-y) > > > if(sum(abs(resgrgap - grcheck)) > tol) + stop("wrong result 2") > > > resgrgap <- gradygapNIR(x0, y, dimx, grobj=grobj) > > grcheck <- c( - grobj(xy1, 1, 1) - grobj(xy2, 2, 1), + - grobj(xy1, 1, 2) - grobj(xy2, 2, 2)) + alpha*(x0-y) > > > if(sum(abs(resgrgap - grcheck)) > tol) + stop("wrong result 3") > > > > > fpNIR(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, echo=FALSE, control=list(eps=1e-10)) $par [1] 0.6610086 0.3389914 $value [1] -0.5337785 $counts function gradient 463 36 $convergence [1] 0 $message NULL $outer.iterations [1] 7 $barrier.value [1] -0.0001723885 $optim.function [1] "constrOptim.nl" $optim.method [1] "BFGS" > fpNIR(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, echo=FALSE, yinit=runif(2, max=1/2)) $par [1] 0.6652873 0.3347127 $value [1] -0.5337201 $counts function gradient 558 29 $convergence [1] 0 $message NULL $outer.iterations [1] 5 $barrier.value [1] -0.00457254 $optim.function [1] "constrOptim.nl" $optim.method [1] "BFGS" > > > > fpNIR(x0, dimx, obj=obj, grobj=grobj, echo=FALSE) $par [1] 0.9634575 0.4882781 $value [1] -0.6790148 $counts function gradient 5 3 $convergence [1] 0 $message NULL $optim.function [1] "optim" $optim.method [1] "BFGS" > fpNIR(x0, dimx, obj=obj, joint=h, echo=FALSE) $par [1] 0.6609286 0.3390714 $value [1] -0.5337785 $counts function gradient 123 27 $convergence [1] 0 $message NULL $outer.iterations [1] 3 $barrier.value [1] -0.0004844439 $optim.function [1] "constrOptim.nl" $optim.method [1] "BFGS" > fpNIR(x0, dimx, obj=obj, echo=FALSE) $par [1] 0.9634416 0.4940028 $value [1] -0.6790818 $counts function gradient 9 6 $convergence [1] 0 $message NULL $optim.function [1] "optim" $optim.method [1] "BFGS" > > res1 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="pure", merit="FP", control.outer=list(maxiter=10, trace=1)) fpiter Iter: 1 Residual: 0.4898871 Iter: 2 Residual: 0.05883914 Iter: 3 Residual: 0.002600994 Iter: 4 Residual: 0.00907679 Iter: 5 Residual: 0.004090908 Iter: 6 Residual: 0.02428721 Iter: 7 Residual: 0.05000278 Iter: 8 Residual: 0.001643868 Iter: 9 Residual: 0.05360592 > res1$inner.counts function gradient 5642 284 > res1$outer.counts fn merit 10 0 > > > res2 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="pure", merit="NI", control.outer=list(maxit=10, echo=2)) **** k 1 x_k 0.6795821 0.3204179 m(x_k) 0.0004777105 **** k 2 x_k 0.6685395 0.3314605 m(x_k) -0.009069673 **** k 3 x_k 0.7208472 0.2791528 m(x_k) 0.004872223 **** k 4 x_k 0.7027293 0.2972707 m(x_k) 0.003766559 **** k 5 x_k 0.6695455 0.3304545 m(x_k) -0.001164071 **** k 6 x_k 0.6863368 0.3136632 m(x_k) -0.002835375 **** k 7 x_k 0.7030362 0.2969638 m(x_k) -0.01573879 **** k 8 x_k 0.7473119 0.2526881 m(x_k) 0.01695409 **** k 9 x_k 0.6942036 0.3057964 m(x_k) 0.0007786858 **** k 10 x_k 0.6889783 0.3110217 m(x_k) 0.0006344455 > res2$inner.counts function gradient 6287 318 > res2$outer.counts fn merit 11 11 > > if(FALSE) + { + res3 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="UR", merit="FP", control.outer=list(maxit=10), stepfunc=decrstep, argstep=5) + + res3 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="UR", merit="FP", control.outer=list(maxit=10, echo=3), stepfunc=decrstep5) + res3$inner.counts + res3$outer.counts + + res4 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="UR", merit="NI", control.outer=list(maxit=10, echo=3), stepfunc=decrstep5) + res4$inner.counts + res4$outer.counts + + res5 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="RRE", merit="NI", control.outer=list(maxit=10, echo=3)) + res5$inner.counts + res5$outer.counts + + res5 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="RRE", merit="FP", control.outer=list(maxiter=10, trace=1), order=1) + res5$inner.counts + res5$outer.counts + + res6 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="MPE", merit="NI", control.outer=list(maxit=10, echo=3)) + res6$inner.counts + res6$outer.counts + + res6 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="MPE", merit="FP", control.outer=list(maxiter=10, trace=1), order=1) + res6$inner.counts + res6$outer.counts + + res7 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="SqMPE", merit="NI", control.outer=list(maxit=10, echo=3)) + res7$inner.counts + res7$outer.counts + + res7 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="SqMPE", merit="FP", control.outer=list(maxiter=10, trace=1), order=1) + res7$inner.counts + res7$outer.counts + + res8 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="SqRRE", merit="NI", control.outer=list(maxit=10, echo=3)) + res8$inner.counts + res8$outer.counts + + res8 <- GNE.fpeq(x0, dimx, obj=obj, grobj=grobj, joint=h, jacjoint=jach, method="SqRRE", merit="FP", control.outer=list(maxiter=10, trace=1)) + res8$inner.counts + res8$outer.counts + } > > > #------------------------------------------------------------------------------- > # (5) another Example > #------------------------------------------------------------------------------- > > # associated objective functions > > dimx <- c(2, 2, 3) > #O_i(x) > fullob <- function(x, i) + { + x <- x[1:7] + if(i == 1) + res <- sum((x - 1:7)^3) + if(i == 2) + res <- sum((x - 1:7)^(1:7)) + if(i == 3) + res <- x[1] + x[3] + (x[5]^2 + x[6]^2 + x[7]^2 - 5)^2 + + res + } > #Gr_x_j O_i(x) > grfullob <- function(x, i, j) + { + x <- x[1:7] + if(i == 1) + grad <- 3*(x - 1:7)^2 + if(i == 2) + grad <- 1:7*(x - 1:7)^(0:6) + if(i == 3) + { + s <- x[5]^2 + x[6]^2 + x[7]^2 - 5 + grad <- c(1, 0, 1, 0, 4*x[5]*s, 4*x[6]*s, 4*x[7]*s) + } + grad[j] + } > > > #Gr_x_k Gr_x_j O_i(x) > hefullob <- function(x, i, j, k) + { + x <- x[1:7] + if(i == 1) + he <- diag( 6*(x - 1:7) ) + if(i == 2) + he <- diag( c(0, 2, 6, 12, 20, 30, 42)*(x - 1:7)^c(0, 0:5) ) + if(i == 3) + { + s <- x[5]^2 + x[6]^2 + x[7]^2 + + he <- rbind(rep(0, 7), rep(0, 7), rep(0, 7), rep(0, 7), + c(0, 0, 0, 0, 4*s+8*x[5]^2, 8*x[5]*x[6], 8*x[5]*x[7]), + c(0, 0, 0, 0, 8*x[5]*x[6], 4*s+8*x[6]^2, 8*x[6]*x[7]), + c(0, 0, 0, 0, 8*x[5]*x[7], 8*x[6]*x[7], 4*s+8*x[7]^2) ) + } + he[j,k] + } > > x0 <- rexp(sum(dimx)) > y <- rexp(sum(dimx)) > xy1 <- c(y[1:2], x0[3:7]) > xy2 <- c(x0[1:2], y[3:4], x0[5:7]) > xy3 <- c(x0[1:4], y[5:7]) > alpha <- 0.1 > > resgap <- gapNIR(x0, y, dimx, obj=fullob) > > gapcheck <- fullob(x0, 1)-fullob(xy1, 1)-alpha/2*sum((x0[1:2]-y[1:2])^2) > gapcheck <- gapcheck + fullob(x0, 2)-fullob(xy2, 2)-alpha/2*sum((x0[3:4]-y[3:4])^2) > gapcheck <- gapcheck + fullob(x0, 3)-fullob(xy3, 3)-alpha/2*sum((x0[5:7]-y[5:7])^2) > > if(sum(abs(resgap - gapcheck)) > tol) + stop("wrong result 4") > > > resgrxgap <- gradxgapNIR(x0, y, dimx, grobj=grfullob) > > grxcheck <- sapply(1:7, function(j) grfullob(x0, 1, j)) - sapply(1:7, function(j) grfullob(xy1, 1, j)) > grxcheck <- grxcheck + sapply(1:7, function(j) grfullob(x0, 2, j)) - sapply(1:7, function(j) grfullob(xy2, 2, j)) > grxcheck <- grxcheck + sapply(1:7, function(j) grfullob(x0, 3, j)) - sapply(1:7, function(j) grfullob(xy3, 3, j)) > grxcheck <- grxcheck + c(sapply(1:2, function(j) grfullob(xy1, 1, j)), sapply(3:4, function(j) grfullob(xy2, 2, j)), sapply(5:7, function(j) grfullob(xy3, 3, j))) > grxcheck <- grxcheck - alpha*(x0-y) > > if(sum(abs(resgrxgap - grxcheck)) > tol) + stop("wrong result 5") > > > resgrygap <- gradygapNIR(x0, y, dimx, grobj=grfullob) > > grycheck <- - c(sapply(1:2, function(j) grfullob(xy1, 1, j)), sapply(3:4, function(j) grfullob(xy2, 2, j)), sapply(5:7, function(j) grfullob(xy3, 3, j))) > grycheck <- grycheck + alpha*(x0-y) > > if(sum(abs(resgrygap - grycheck)) > tol) + stop("wrong result 6") > > > > proc.time() user system elapsed 3.42 0.15 3.59