spg Projection Functions

Description

Projection function implementing contraints for spg parameters.

Usage

     projectLinear(par, A, b, meq)

Arguments

Details

The function projectLinear can be used by spg to define the constraints of the problem. It projects a point in R^n onto a region that defines the constraints. It takes a real vector par as argument and returns a real vector of the same length.

The function projectLinear incorporates linear equalities and inequalities in nonlinear optimization using a projection method, where an infeasible point is projected onto the feasible region using a quadratic programming solver. The inequalities are defined such that: A %*% x - b > 0 . The first ‘meq’ rows of A and the first ‘meq’ elements of b correspond to equality constraints.

Value

A vector of the constrained parameter values.

Note

We are grateful to Berwin Turlach for creating a special version of solve.QP function that is specifically tailored to solving the projection problem.

See Also

spg

Examples

# Example
fn <- function(x) (x[1] - 3/2)^2 + (x[2] - 1/8)^4

gr <- function(x) c(2 * (x[1] - 3/2) , 4 * (x[2] - 1/8)^3)

# This is the set of inequalities
# x[1] - x[2] >= -1
# x[1] + x[2] >= -1
# x[1] - x[2] <= 1
# x[1] + x[2] <= 1

# The inequalities are written in R such that:  Amat %*% x  >= b 
Amat <- matrix(c(1, -1, 1, 1, -1, 1, -1, -1), 4, 2, byrow=TRUE)
b <- c(-1, -1, -1, -1)
meq <- 0  # all 4 conditions are inequalities

p0 <- rnorm(2)
spg(par=p0, fn=fn, gr=gr, project="projectLinear", 
      projectArgs=list(A=Amat, b=b, meq=meq))

meq <- 1  # first condition is now an equality
spg(par=p0, fn=fn, gr=gr, project="projectLinear", 
      projectArgs=list(A=Amat, b=b, meq=meq))


# box-constraints can be incorporated as follows:
# x[1] >= 0
# x[2] >= 0
# x[1] <= 0.5
# x[2] <= 0.5

Amat <- matrix(c(1, 0, 0, 1, -1, 0, 0, -1), 4, 2, byrow=TRUE)
b <- c(0, 0, -0.5, -0.5)

meq <- 0
spg(par=p0, fn=fn, gr=gr, project="projectLinear", 
   projectArgs=list(A=Amat, b=b, meq=meq))

# Note that the above is the same as the following:
spg(par=p0, fn=fn, gr=gr, lower=0, upper=0.5)


# An example showing how to impose other constraints in spg()

fr <- function(x) { ## Rosenbrock Banana function
  x1 <- x[1] 
  x2 <- x[2] 
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2 
  } 

# Impose a constraint that sum(x) = 1

proj <- function(x){ x / sum(x) }

spg(par=runif(2), fn=fr, project="proj") 

# Illustration of the importance of `projecting' the constraints, rather 
#   than simply finding a feasible point:

fr <- function(x) { ## Rosenbrock Banana function 
x1 <- x[1] 
x2 <- x[2] 
100 * (x2 - x1 * x1)^2 + (1 - x1)^2 
} 
# Impose a constraint that sum(x) = 1 

proj <- function(x){ 
# Although this function does give a feasible point it is 
#  not a "projection" in the sense of the nearest feasible point to `x'
x / sum(x) 
} 

p0 <- c(0.93, 0.94)  

# Note, the starting value is infeasible so the next 
#   result is "Maximum function evals exceeded"

spg(par=p0, fn=fr, project="proj") 

# Correct approach to doing the projection using the `projectLinear' function

spg(par=p0, fn=fr, project="projectLinear", projectArgs=list(A=matrix(1, 1, 2), b=1, meq=1)) 

# Impose additional box constraint on first parameter

p0 <- c(0.4, 0.94)    # need feasible starting point

spg(par=p0, fn=fr,  lower=c(-0.5, -Inf), upper=c(0.5, Inf),
  project="projectLinear", projectArgs=list(A=matrix(1, 1, 2), b=1, meq=1)) 

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

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.

> library(BB)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/BB/project.Rd_%03d_medium.png", width=480, height=480)
> ### Name: project
> ### Title: spg Projection Functions
> ### Aliases: projectLinear
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> # Example
> fn <- function(x) (x[1] - 3/2)^2 + (x[2] - 1/8)^4
> 
> gr <- function(x) c(2 * (x[1] - 3/2) , 4 * (x[2] - 1/8)^3)
> 
> # This is the set of inequalities
> # x[1] - x[2] >= -1
> # x[1] + x[2] >= -1
> # x[1] - x[2] <= 1
> # x[1] + x[2] <= 1
> 
> # The inequalities are written in R such that:  Amat %*% x  >= b 
> Amat <- matrix(c(1, -1, 1, 1, -1, 1, -1, -1), 4, 2, byrow=TRUE)
> b <- c(-1, -1, -1, -1)
> meq <- 0  # all 4 conditions are inequalities
> 
> p0 <- rnorm(2)
> spg(par=p0, fn=fn, gr=gr, project="projectLinear", 
+       projectArgs=list(A=Amat, b=b, meq=meq))
iter:  0  f-value:  11.99853  pgrad:  1.859255 
$par
[1] 1 0

$value
[1] 0.2502441

$gradient
[1] 5.551115e-17

$fn.reduction
[1] 11.74829

$iter
[1] 1

$feval
[1] 2

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> meq <- 1  # first condition is now an equality
> spg(par=p0, fn=fn, gr=gr, project="projectLinear", 
+       projectArgs=list(A=Amat, b=b, meq=meq))
iter:  0  f-value:  11.99853  pgrad:  0.8592554 
$par
[1] 0 1

$value
[1] 2.836182

$gradient
[1] 0

$fn.reduction
[1] 9.162349

$iter
[1] 1

$feval
[1] 2

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> 
> # box-constraints can be incorporated as follows:
> # x[1] >= 0
> # x[2] >= 0
> # x[1] <= 0.5
> # x[2] <= 0.5
> 
> Amat <- matrix(c(1, 0, 0, 1, -1, 0, 0, -1), 4, 2, byrow=TRUE)
> b <- c(0, 0, -0.5, -0.5)
> 
> meq <- 0
> spg(par=p0, fn=fn, gr=gr, project="projectLinear", 
+    projectArgs=list(A=Amat, b=b, meq=meq))
iter:  0  f-value:  11.99853  pgrad:  0.5 
$par
[1] 0.5000000 0.1145664

$value
[1] 1

$gradient
[1] 4.543266e-06

$fn.reduction
[1] 10.99853

$iter
[1] 7

$feval
[1] 8

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> # Note that the above is the same as the following:
> spg(par=p0, fn=fn, gr=gr, lower=0, upper=0.5)
iter:  0  f-value:  11.99853  pgrad:  0.5 
$par
[1] 0.5000000 0.1145664

$value
[1] 1

$gradient
[1] 4.543266e-06

$fn.reduction
[1] 10.99853

$iter
[1] 7

$feval
[1] 8

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> 
> # An example showing how to impose other constraints in spg()
> 
> fr <- function(x) { ## Rosenbrock Banana function
+   x1 <- x[1] 
+   x2 <- x[2] 
+   100 * (x2 - x1 * x1)^2 + (1 - x1)^2 
+   } 
> 
> # Impose a constraint that sum(x) = 1
> 
> proj <- function(x){ x / sum(x) }
> 
> spg(par=runif(2), fn=fr, project="proj") 
iter:  0  f-value:  9.888397  pgrad:  0.2925466 
$par
[1] 0.6186303 0.3813697

$value
[1] 0.1456207

$gradient
[1] 1.533979e-06

$fn.reduction
[1] 9.742776

$iter
[1] 8

$feval
[1] 9

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> # Illustration of the importance of `projecting' the constraints, rather 
> #   than simply finding a feasible point:
> 
> fr <- function(x) { ## Rosenbrock Banana function 
+ x1 <- x[1] 
+ x2 <- x[2] 
+ 100 * (x2 - x1 * x1)^2 + (1 - x1)^2 
+ } 
> # Impose a constraint that sum(x) = 1 
> 
> proj <- function(x){ 
+ # Although this function does give a feasible point it is 
+ #  not a "projection" in the sense of the nearest feasible point to `x'
+ x / sum(x) 
+ } 
> 
> p0 <- c(0.93, 0.94)  
> 
> # Note, the starting value is infeasible so the next 
> #   result is "Maximum function evals exceeded"
> 
> spg(par=p0, fn=fr, project="proj") 
iter:  0  f-value:  0.568901  pgrad:  0.002673382 
$par
[1] 0.4973262 0.5026738

$value
[1] 0.568901

$gradient
[1] 0.003780733

$fn.reduction
[1] 0

$iter
[1] 1

$feval
[1] 1

$convergence
[1] 2

$message
[1] "Maximum function evals exceeded"

Warning message:
In spg(par = p0, fn = fr, project = "proj") : Unsuccessful convergence.
> 
> # Correct approach to doing the projection using the `projectLinear' function
> 
> spg(par=p0, fn=fr, project="projectLinear", projectArgs=list(A=matrix(1, 1, 2), b=1, meq=1)) 
iter:  0  f-value:  0.568901  pgrad:  52.24003 
$par
[1] 0.6187956 0.3812044

$value
[1] 0.145607

$gradient
[1] 1.670886e-07

$fn.reduction
[1] 0.423294

$iter
[1] 5

$feval
[1] 9

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> # Impose additional box constraint on first parameter
> 
> p0 <- c(0.4, 0.94)    # need feasible starting point
> 
> spg(par=p0, fn=fr,  lower=c(-0.5, -Inf), upper=c(0.5, Inf),
+   project="projectLinear", projectArgs=list(A=matrix(1, 1, 2), b=1, meq=1)) 
iter:  0  f-value:  61.2  pgrad:  0.27 
$par
[1] 0.5000000 0.4999999

$value
[1] 6.499997

$gradient
[1] 5.960464e-08

$fn.reduction
[1] 54.7

$iter
[1] 1

$feval
[1] 2

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>