Introduction

This vignette goes beyond the Introduction to the levelSets Package vignette to provide more details about how to set up and run a level set problem. It uses a simple 2D problem for illustration. More extended examples are provided in the Example vignettes.

Recall that a level set of a response function \(f(x)\) is the set of \(d\)-dimensional input points \(x\) for which the function value is greater than or equal to a specified threshold \(t\): \[L_f(t) = \{x: f(x) \ge t\}\] If there are restrictions on the inputs that \(f\) will accept, we call the allowed portion of the input space the feasible region for \(f\), and any point outside the feasible region is infeasible.

This package maps out the boundary of a level set by finding its intersections with collections of 1-dimensional rays. A ray is just a line with a defined origin (referred to as position 0 of the ray) and orientation. A ray’s orientation is specified by a second, distinct point that defines position 1 along the ray. Any other point on the ray has a position value obtained by interpolation or extrapolation from the origin and position 1 points. (Ray points on the opposite side of the origin from the position 1 point have negative position values.)

The package provides tools to generate rays, find intersections, and visualize results.

An example

The Rosenbrock “banana” function will be used for illustration. It has a 2-dimensional input space, and a minimum value of 0 at the point (1, 1). Since level sets are defined to have function values at or above a threshold, we use the negative of the function: \[f(x_1, x_2) = -100((x_2 - x_1^2)^2) - (1 - x_1)^2\] Contours of this function can be strongly curved, so level sets may not be convex. The level set threshold we use is \(f \ge -50\).

The response function can be coded in R as

respf <- function(x, inclGrad=FALSE) {
  y <- -1 * (100*((x[, 2] - x[, 1]^2)^2) + (1 - x[, 1])^2)
  if (inclGrad) {
    grad <- cbind(400*(x[, 1]*(x[, 2] - x[, 1]^2)) + 2*(1 - x[, 1]), 
                  -200*(x[, 2] - x[, 1]^2))
    attr(y, "gradient") <- grad
  }
  y
}
thresh <- -50

The gradient calculation is simple so it is included, but gradients are not necessary in general.

Although this function is well-defined for any \((x_1, x_2)\), for illustration we define the feasible region to be \(\{x: x_1 \ge -2\}\). A feasibility function can be coded as

feasf <- function(x) {
     (x[, 1] >= -2)
}

In this package point coordinates are always represented as rows of a \(d\)-column matrix, one row per point. Both the response and feasibility functions must accept as their first argument a matrix of point coordinates, and return a vector with one value per point: a numeric value for the response function, and TRUE or FALSE for the feasibility function.

The following sections describe package functions and classes, in the order they are likely to be used in studying a level set.

Specifying the response function and its input space

library(levelSets)
fobj <- fnObj(c("x1", "x2"), respfn=respf, feasfn=feasf, hasgrad=TRUE)

str(lsetPkgOpt(c("inptol", "resptol")))

## List of 2
##  $ inptol : num [1:2] 1.00e-06 1.49e-08
##  $ resptol: num [1:2] 1.00e-05 1.49e-08

Specifying the search region for a level set

The user must specify a finite region of input space where the search for level set boundary points will be focused. The search region takes the form of a \(d\)-dimensional box or hyper-rectangle, represented by objects of class inpRect, and created by the inpRect() function. See the documentation for the latter for details.

Ideally the search region should be slightly larger than the smallest hyper-rectangle that contains the level set. However if an initial guess at the search region turns out to be too small, too large, or simply in the wrong part of input space, the search can be resumed with an updated search region. Functions boundingRect() and segBoundingRect() may help with choosing that.

Note that the search region is distinct from the feasible region of the response function. The search region is allowed to include infeasible points. In fact, when the level set touches the boundary of the feasible region, it actually helps the search algorithm if the search region extends beyond the feasible region in that direction.

An example search region for the banana function:

rect0 <- inpRect(cbind(c(-3, 3), c(-3, 6)), spec=fobj)

Specifying search rays

Points on the boundary of a level set are identified by finding where 1-dimensional rays intersect the boundary. Collections of rays, with a shared origin, are represented by inpRays objects. As mentioned in the Introduction, points along a ray are identified by their position relative to the ray origin (position 0) and second defining point for the ray (position 1). Rays are in principle infinite so there are points at any position value, positive or negative. There is no requirement that the Euclidean distance between positions 0 and 1 be the same for different rays in an inpRays object.

Degenerate rays, where the position 0 and position 1 points are equal or nearly equal, are not allowed. Whether two points are considered nearly equal is determined by the numerical tolerances for point coordinates, as described for inptol in section 3.

A set of rays can be created with the inpRays() function. The randomRays() and axisRays() functions create rays with specific properties; see their documentation.

set.seed(1)
rays0 <- randomRays(25, origin=c(1, 1), rect=rect0)
plot(rays0, main="Random rays, position 1 defined by 'rect0'")

Scaling of input space dimensions

For some problems the dimensions of the input space have different units, or point coordinates for different dimensions have different magnitudes, or the level set of interest is elongated in certain directions. In these cases it is important to specify that package algorithms should rescale and/or weight the dimensions when they calculate distances or angles in input space. inpScale objects provide a way to specify such scaling. (This is different from the redefinition of position 1 along rays discussed in the previous section.)

The following example in two dimensions illustrates the issue. Suppose the level set is believed to have a wider range of \(x_1\) values than \(x_2\) values, and the search region reflects this:

fspec <- fnSpec(inpnames=c("x1", "x2"))
rect1 <- inpRect(rbind(c(0.5, 3), c(5.5, 3.5)), spec=fspec)

By default randomRays() creates ray vectors that are uniformly distributed around the unit circle; i.e. it treats both dimensions equally and symmetrically. When such rays are used with an asymmetric search region we can end up undersampling the long dimension of the region, as seen in the first plot below. Scaling each dimension in inverse proportion to the lengths of the sides of the search region gives a better distribution (second plot).

set.seed(1)
# Rays with uniformly distributed directions and equal position 1 lengths:
rays_u <- randomRays(50, origin=c(1, 3.25), spec=fspec)
# Position 1 redefined by intersection with a hyper-rectangle:
rays1 <- update(rays_u, rect=rect1)
plot(rays1, main="Random rays, no dimension scaling")
# Specify scaling for input space dimensions:
scl2 <- inpScale(rectSize(rect1), spec=fspec)
# Same rectangle, but ray generation now reflects input dimension scaling:
set.seed(1)
rays2 <- randomRays(50, origin=c(1, 3.25), scale=scl2, rect=rect1)
plot(rays2, main="Random rays, dimensions scaled by 'rect' sides")

Scaling factors should be chosen so that division by them standardizes point coordinates: coordinates for different dimensions should have roughly similar magnitude and spread, preferably with a magnitude not too far from 1.

inpScale objects in fact allow for more general types of scaling, such as correlation between input space dimensions. See ?inpScale for details.

Profiling a response function along rays

The respProfiles() function allows one to evaluate the response function at an arbitrary set of positions along each of a set of rays. Returning to the banana response function of section 2:

rays3 <- axisRays(degree=2, origin=c(1, 1), rect=rect0)
prof3 <- respProfiles(fobj, rays=rays3, lsetthresh=-50, posns=(-10:10)/10)

axisRays() with degree=2 generates four rays: one parallel to each of the two coordinate axes, and two parallel to the diagonals with slopes of \(\pm 1\). The posns argument to respProfiles specifies the evaluation positions (the same positions are used for each ray). The optional lsetthresh argument specifies the level set of interest. It triggers a search for level set boundary points along each ray within the range of posns. There are summary and plot methods for the returned respProfiles object. (The plot method requires the ggplot2 package.)

summary(prof3)

## Response function profiled along 4 rays in a 2-dimensional input space
## Level set threshold: -50
## Number of profile points over all rays: 93
##   In feasible region= 93 (on its bdry= 0)
##   In level set= 32 (on its bdry= 9)
## Range of response values over all rays: -10004, 0

if (requireNamespace("ggplot2", quietly=TRUE))  plot(prof3, groupBy=names(rays3))

Function respInfo() applied to a respProfiles object returns a data frame with information about response at each of the positions along each ray. Function ptCoord() returns a matrix with the corresponding point coordinates.

head(respInfo(prof3))

##        resp feas  lset line       posn feasbdry lsetbdry    deriv      trend
## 1  -4.00000 TRUE  TRUE    1 -1.0000000       NA       NA    8.000 increasing
## 2 -16.20000 TRUE  TRUE    1 -0.9000000    FALSE    FALSE -223.200 decreasing
## 3 -43.52000 TRUE  TRUE    1 -0.8000000    FALSE    FALSE -300.800 decreasing
## 4 -49.99996 TRUE  TRUE    1 -0.7784998    FALSE     TRUE -301.125 decreasing
## 5 -72.52000 TRUE FALSE    1 -0.7000000    FALSE    FALSE -263.200 decreasing
## 6 -93.60000 TRUE FALSE    1 -0.6000000    FALSE    FALSE -148.800 decreasing

head(ptCoord(prof3))

##              x1 x2
## [1,] -1.0000000  1
## [2,] -0.8000000  1
## [3,] -0.6000000  1
## [4,] -0.5569996  1
## [5,] -0.4000000  1
## [6,] -0.2000000  1

Level set boundary search

The emphasis of respProfiles() is evaluation of the response function at a specified set of ray positions. Finding level set boundary points is optional and secondary. When finding boundary points is the primary goal, other functions are available.

slices <- rbind(c(1, NA, NA), c(2, NA, 3), c(NA, NA, NA))