project-onto-variety: Projection onto a variety

In dkahle/algstat: Algebraic Statistics

Projection onto a variety

Description

A R-based implementation of the gradient descent homotopies. The lagrange version uses Newton's method on the Lagrangian system.

Usage

project_onto_variety(
  x0,
  poly,
  dt = 0.05,
  varorder = sort(vars(poly)),
  n_correct = 2,
  al = rnorm(length(x0)),
  message = FALSE,
  tol = .Machine$double.eps^(1/2),
  gfunc,
  dgfunc,
  ddgfunc,
  bias = 0
)

project_onto_variety_lagrange(
  x0,
  poly,
  varorder = vars(poly),
  method = "newton",
  maxit = 1000,
  tol = .Machine$double.eps^(1/2),
  tol_x = .Machine$double.eps^(1/2),
  message = FALSE,
  ...
)

project_onto_variety_gradient_descent(
  x0,
  poly,
  varorder = vars(poly),
  ga = 0.01,
  max_ga = 0.1,
  method = c("line", "optimal", "fixed"),
  tol = .Machine$double.eps^(1/2),
  tol_x = .Machine$double.eps^(1/2),
  maxit = 1000,
  message = FALSE
)

project_onto_variety_newton(
  x0,
  poly,
  varorder = vars(poly),
  ga = 1e-04,
  max_ga = 2,
  method = c("line", "fixed"),
  tol = .Machine$double.eps^(1/2),
  tol_x = .Machine$double.eps^(1/2),
  maxit = 1000,
  message = FALSE
)

Arguments

x0	Atomic vector, the point to be projected.
poly	An mpoly object, typically created with mp().
dt	The t-mesh size for the homotopy.
varorder	A character vector specifying the variable order to pass to mpoly::as.function.mpoly().
n_correct	The number of Newton correction iterations to use.
al	A numeric vector of length 2; the patch to do projective calculations over.
message	If TRUE, the user is issued messages on the algorithm's progress.
tol	A tolerance on the residual; a warning is issued if the magnitude of the residual is larger than tol.
gfunc, dgfunc, ddgfunc	The polynomial poly, its gradient, and its Hessian as functions. Only used in project_onto_variety(), and computed internally if not provided.
bias	A multiple to add to the identity to make the Jacobian invertible.
method	Used in project_onto_variety_lagrange() and project_onto_variety_gradient_descent(). In the former, if "newton", a simple R implementation of Newton's method to solve the nonlinear algebraic system generated by the Lagrangian; otherwise, a character string to pass to optim() to minimize the sum of the squares of the Lagrangian. In the latter, the method of selecting the learning rate, "line" (line search) or "fixed".
maxit	Number of maximum iterations in solving Newton's method.
tol_x	A tolerance on subsequent step sizes.
...	Additional arguments to pass to optim() when method is not "newton".
ga	Learning rate for gradient descent.
max_ga	Maximum learning rate for gradient descent line search.

Value

A numeric vector the same length as x0.

Author(s)

David Kahle

References

Griffin, Z. and J. Hauenstein (2015). Real solutions to systems of polynomial equations and parameter continuation. Advances in Geometry 15(2), pp.173–187.

Bates, D., J. Hauenstein, A. Sommese, and C. Wampler (2013). Numerically Solving Polynomial Systems with Bertini. SIAM. pp.34–35.

Examples

library("ggplot2")


## basic usage
########################################

x0 <- c(1,1)
p <- mp("x^2 + y^2 - 1")
(x0_proj <- project_onto_variety(x0, p))

as.function(p)(x0_proj)
sqrt(2)/2

cbind(t(x0), t(x0_proj)) %>%
  as.data.frame() %>% tibble::as_tibble() %>%
  purrr::set_names(c("x", "y", "x_proj", "y_proj")) -> df

ggvariety(p) + coord_equal() +
  geom_segment(
    aes(x, y, xend = x_proj, yend = y_proj),
    data = df, inherit.aes = FALSE
  )


# alternatives
1 / sqrt(2)
project_onto_variety_lagrange(x0, p)
project_onto_variety_newton(x0, p)
project_onto_variety_gradient_descent(x0, p)
project_onto_variety_gradient_descent(x0, p, method = "line")


# number of variables > 2
x0 <- c(1,1,1)
p <- mp("x^2 + y^2 + z^2 - 1")
1 / sqrt(3)
project_onto_variety(x0, p)
project_onto_variety_lagrange(x0, p)
project_onto_variety_newton(x0, p)
project_onto_variety_gradient_descent(x0, p)
project_onto_variety_gradient_descent(x0, p, method = "line")


## options
########################################

x0 <- c(1,1)
p <- mp("x^2 + y^2 - 1")
project_onto_variety(x0, p, message = TRUE)
project_onto_variety(x0, p, dt = .25, message = TRUE)


# precomputing the function, gradient, and hessian
varorder <- c("x", "y")
gfunc <- as.function(p, varorder = varorder)

dg <- deriv(p, var = varorder)
dgfunc <- as.function(dg, varorder = varorder)

ddg <- lapply(dg, deriv, var = varorder)
ddgfunc_list <- lapply(ddg, as.function, varorder = varorder, silent = TRUE)
ddgfunc <- function(x) sapply(ddgfunc_list, function(f) f(x))

project_onto_variety(x0, p, gfunc = gfunc, dgfunc = dgfunc, ddgfunc = ddgfunc)


## more complex example
########################################

p <- mp("(x^2 + y^2)^2 - 2 (x^2 - y^2)")
ggvariety(p, c(-2, 2), n = 201) + coord_equal()

x0 <- c(.025, .30)
(x0_proj <- project_onto_variety(x0, p))

cbind(t(x0), t(x0_proj)) %>%
  as.data.frame() %>% tibble::as_tibble() %>%
  purrr::set_names(c("x", "y", "x_proj", "y_proj")) -> df

ggvariety(p, c(-2, 2)) + coord_equal() +
  geom_segment(
    aes(x, y, xend = x_proj, yend = y_proj),
    data = df, inherit.aes = FALSE
  )




## projecting a dataset - grid
########################################

library("ggplot2")
library("dplyr")

(p <- lissajous(5, 5, 0, 0))
# (p <- lissajous(9, 9, 0, 0))
# p <- mp("x^2 + y^2 - 1")
ggvariety(p, n = 251) + coord_equal()

set.seed(1)
(s <- seq(-1, 1, .25))
n <- length(s)
grid <- expand.grid(x = s, y = s)
grid$x <- jitter(grid$x)
grid$y <- jitter(grid$y)

ggplot(grid, aes(x, y)) + geom_point() + coord_equal()

grid_proj <- project_onto_variety(grid, p)
head(grid_proj)
names(grid_proj) <- c("x_proj", "y_proj")

ggvariety(p, n = 251) + coord_equal() +
  geom_segment(
    aes(x, y, xend = x_proj, yend = y_proj),
    data = bind_cols(grid, grid_proj), inherit.aes = FALSE
  ) +
  geom_point(aes(x, y), data = grid, inherit.aes = FALSE)


# here's what happens when you use a naive implementation -
# gradient descent on g^2 with line search
grid_proj_gd <- project_onto_variety_gradient_descent(grid, p, method = "optimal")
names(grid_proj_gd) <- c("x_proj", "y_proj")

grid_proj_lagrange <- project_onto_variety_lagrange(grid, p)
names(grid_proj_lagrange) <- c("x_proj", "y_proj")

grid_proj_newton <- project_onto_variety_newton(grid, p)
names(grid_proj_newton) <- c("x_proj", "y_proj")

df <- bind_rows(
  bind_cols(grid, grid_proj) %>% mutate(method = "gradient descent homotopy"),
  bind_cols(grid, grid_proj_gd) %>% mutate(method = "optimal gradient descent"),
  bind_cols(grid, grid_proj_newton) %>% mutate(method = "newton"),
  bind_cols(grid, grid_proj_lagrange) %>% mutate(method = "newton on lagrangian")
)

ggvariety(p, n = 251) +
  geom_segment(
    aes(x, y, xend = x_proj, yend = y_proj),
    data = df, inherit.aes = FALSE
  ) +
  geom_point(aes(x, y), data = grid, inherit.aes = FALSE) +
  coord_equal() +
  facet_wrap(~ method) + xlim(-1.1, 1.1) + ylim(-1.1, 1.1)


## projecting a dataset - rvnorm
########################################

library("ggplot2")
library("dplyr")

## Not run:  requires stan

(p <- lissajous(5, 5, 0, 0))
ggvariety(p, n = 251) + coord_equal()

set.seed(1)
(samps <- rvnorm(1e4, p, sd = .025, output = "tibble"))

ggplot(samps, aes(x, y)) +
  geom_point(aes(color = chain)) +
  coord_equal() +
  facet_wrap(~ chain)

ggplot(samps, aes(x, y)) +
  geom_bin2d(binwidth = .03*c(1,1)) +
  coord_equal()

# cut down on draws for time
subsamps <- samps %>% sample_n(500)
ggplot(subsamps, aes(x, y)) + geom_point() + coord_equal()

subsamps %>%
  select(x, y) %>%
  as.matrix() %>%
  #apply(1, function(x0) project_onto_variety(x0, p, dt = .025, n_correct = 3)) %>% t() %>%
  apply(1, function(x0) project_onto_variety_lagrange(x0, p)) %>% t() %>%
  as.data.frame() %>% tibble::as_tibble() %>%
  purrr::set_names(c("x_proj", "y_proj")) %>%
  bind_cols(subsamps, .) ->
  subsamps

ggvariety(p, n = 251) + coord_equal() +
  geom_segment(
    aes(x, y, xend = x_proj, yend = y_proj),
    data = subsamps, inherit.aes = FALSE
  ) +
  geom_point(
    aes(x, y, color = factor(chain)),
    data = subsamps, inherit.aes = FALSE
  )

ggplot(subsamps, aes(x_proj, y_proj)) + geom_point() + coord_equal()


## End(Not run)

dkahle/algstat documentation built on May 23, 2023, 12:29 a.m.