tests/testthat/test-g01-grad.R In anqif/cvxr: Disciplined Convex Optimization

```context("test-g01-grad")
a <- Variable(name = "a")

x <- Variable(2, name = "x")
y <- Variable(2, name = "y")

A <- Variable(2, 2, name = "A")
B <- Variable(2, 2, name = "B")
C <- Variable(3, 2, name = "C")

value(C) <- rbind(c(1,-2), c(3,4), c(-1,-3))
value(A) <- rbind(c(3,2), c(-5,1))
expect_warning(expr <- C %*% A)

c(-5,0,0,1,0,0), c(0,-5,0,0,1,0), c(0,0,-5,0,0,1)))
c(0,0,0,-2,4,-3)))
})

value(x) <- c(-1,0)
expr <- p_norm(x, 1)

value(x) <- c(0,10)
expr <- p_norm(x, 1)

value(x) <- c(-3,4)
expr <- p_norm(x, 2)

value(x) <- c(-1,2)
expr <- p_norm(x, 0.5)

value(x) <- c(0,0)
expr <- p_norm(x, 0.5)

value(x) <- c(0,0)
expr <- p_norm(x, 2)

value(A) <- rbind(c(2,-2), c(2,2))
expr <- p_norm(A, 2)

value(A) <- rbind(c(3,-3), c(4,4))
expr <- p_norm(A, 2, axis = 2)

value(A) <- rbind(c(3,-4), c(4,3))
expr <- p_norm(A, 2, axis = 1)

value(A) <- rbind(c(3,-4), c(4,3))
expr <- p_norm(A, 0.5)
})

value(x) <- c(0,1)
expr <- log_sum_exp(x)
e <- exp(1)

value(A) <- rbind(c(0,1), c(-1,0))
expr <- log_sum_exp(A)

value(A) <- rbind(c(0,1), c(-1,0))
expr <- log_sum_exp(A, axis = 2)
expect_equivalent(as.matrix(grad(expr)[[as.character(A@id)]]), cbind(c(1.0/(1+1.0/e), 1.0/e/(1+1.0/e), 0, 0), c(0, 0, e/(1+e), 1.0/(1+e))))
})

value(x) <- c(1,2)
expr <- geo_mean(x)

value(x) <- c(0,2)
expr <- geo_mean(x)

value(x) <- c(1,2)
expr <- geo_mean(x, c(1,0))

# No exception for single weight
value(x) <- c(-1,2)
expr <- geo_mean(x, c(1,0))
})

value(A) <- cbind(c(2,0), c(0,1))
expr <- lambda_max(A)

value(A) <- cbind(c(1,0), c(0,2))
expr <- lambda_max(A)

value(A) <- cbind(c(1,0), c(0,1))
expr <- lambda_max(A)
})

value(A) <- diag(2)
value(B) <- diag(2)
expr <- matrix_frac(A, B)

value(B) <- matrix(0, nrow = 2, ncol = 2)
expr <- matrix_frac(A, B)

value(x) <- c(2,3)
value(A) <- diag(2)
expr <- matrix_frac(x, A)
})

value(A) <- cbind(c(10,4), c(4,30))
expr <- norm_nuc(A)
})

value(A) <- 2*diag(rep(1,2))
expr <- log_det(A)

mat <- rbind(c(1,2), c(3,5))
value(A) <- t(mat) %*% mat
expr <- log_det(A)
val <- t(base::solve(value(A)))

value(A) <- matrix(0, nrow = 2, ncol = 2)
expr <- log_det(A)

value(A) <- -rbind(c(1,2), c(3,4))
expr <- log_det(A)
})

value(x) <- c(1,2)
value(a) <- 2

value(a) <- 0

value(A) <- diag(rep(1,2))
value(a) <- 2

value(x) <- c(1,2)
value(a) <- 2
value(y) <- c(1,2)
value(a) <- 2
})

value(x) <- c(2,1)
expr <- max_entries(x)

value(A) <- rbind(c(1,2), c(4,3))
expr <- max_entries(A)

value(A) <- rbind(c(1,2), c(4,3))
expr <- max_entries(A, axis = 2)

value(A) <- rbind(c(1,2), c(4,3))
expr <- max_entries(A, axis = 1)
})

test_that("Test sigma_max", {
value(A) <- cbind(c(1,0), c(0,2))
expr <- sigma_max(A)

value(A) <- cbind(c(1,0), c(0,1))
expr <- sigma_max(A)
})

test_that("Test sum_largest", {
value(A) <- cbind(c(4,3), c(2,1))
expr <- sum_largest(A, 2)

value(A) <- cbind(c(1,2), c(3,0.5))
expr <- sum_largest(A, 2)
})

test_that("Test abs", {
value(A) <- cbind(c(1,2), c(-1,0))
expr <- abs(A)
val <- diag(c(1,1,-1,0))
})

test_that("Test linearize method", {
# Affine
value(x) <- c(1,2)
expr <- (2*x - 5)[1]
lin_expr <- linearize(expr)
value(x) <- c(55,22)
expr <- (2*x - 5)[1]
lin_expr <- linearize(expr)
expect_equal(value(lin_expr), value(expr))
value(x) <- c(-1,-5)
expr <- (2*x - 5)[1]
lin_expr <- linearize(expr)
expect_equal(value(lin_expr), value(expr))

# Convex
expr <- A^2 + 5
expect_error(linearize(expr))

value(A) <- cbind(c(1,2), c(3,4))
expr <- A^2 + 5
lin_expr <- linearize(expr)
manual <- value(expr) + 2*reshape_expr(value(diag(vec(A))) %*% vec(A - value(A)), c(2, 2))
expect_equivalent(as.matrix(value(lin_expr)), value(expr))

value(A) <- rbind(c(-5,-5), c(8.2,4.4))
expr <- A^2 + 5
lin_expr <- linearize(expr)
manual <- value(expr) + 2*reshape_expr(value(diag(vec(A))) %*% vec(A - value(A)), c(2, 2))
expect_true(all(value(lin_expr) <= value(expr)))
expect_equivalent(as.matrix(value(lin_expr)), value(manual))

# Concave
value(x) <- c(1,2)
expr <- log(x)/2
lin_expr <- linearize(expr)
manual <- value(expr) + value(diag(0.5*x^-1)) %*% (x - value(x))
expect_equivalent(as.matrix(value(lin_expr)), value(expr))

value(x) <- c(3,4.4)
expr <- log(x)/2
lin_expr <- linearize(expr)
manual <- value(expr) + value(diag(0.5*x^-1)) %*% (x - value(x))
expect_true(all(value(lin_expr) >= value(expr)))
expect_equivalent(as.matrix(value(lin_expr)), value(manual))
})

value(a) <- 2
expr <- log(a)

value(a) <- 3
expr <- log(a)

value(a) <- -1
expr <- log(a)

value(x) <- c(3,4)
expr <- log(x)
val <- diag(c(1/3,1/4))

value(x) <- c(-1e-9,4)
expr <- log(x)

value(A) <- cbind(c(1,2), c(3,4))
expr <- log(A)
val <- diag(c(1, 1/2, 1/3, 1/4))
})

test_that("Test domain for log1p", {
value(a) <- 2
expr <- log1p(a)

value(a) <- 3
expr <- log1p(a)

value(a) <- -1
expr <- log1p(a)

value(x) <- c(3,4)
expr <- log1p(x)
val <- diag(c(1/4,1/5))

value(x) <- c(-1e-9-1,4)
expr <- log1p(x)

value(A) <- cbind(c(1,2), c(3,4))
expr <- log1p(A)
val <- diag(c(1/2, 1/3, 1/4, 1/5))
})

test_that("Test domain for entr", {
value(a) <- 2
expr <- entr(a)

value(a) <- 3
expr <- entr(a)

value(a) <- -1
expr <- entr(a)

value(x) <- c(3,4)
expr <- entr(x)
val <- diag(-(log(c(3,4)) + 1))

value(x) <- c(-1e-9,4)
expr <- entr(x)

value(A) <- cbind(c(1,2), c(3,4))
expr <- entr(A)
val <- diag(-(log(1:4)+1))
})

test_that("Test domain for exp", {
value(a) <- 2
expr <- exp(a)

value(a) <- 3
expr <- exp(a)

value(a) <- -1
expr <- exp(a)

value(x) <- c(3,4)
expr <- exp(x)
val <- diag(exp(c(3,4)))

value(x) <- c(-1e-9,4)
expr <- exp(x)
val <- diag(exp(c(-1e-9,4)))

value(A) <- cbind(c(1,2), c(3,4))
expr <- exp(A)
val <- diag(exp(1:4))
})

test_that("Test domain for logistic", {
value(a) <- 2
expr <- logistic(a)

value(a) <- 3
expr <- logistic(a)

value(a) <- -1
expr <- logistic(a)

value(x) <- c(3,4)
expr <- logistic(x)
val <- diag(exp(c(3,4))/(1+exp(c(3,4))))

value(x) <- c(-1e-9,4)
expr <- logistic(x)
val <- diag(exp(c(-1e-9,4))/(1+exp(c(-1e-9,4))))

value(A) <- cbind(c(1,2), c(3,4))
expr <- logistic(A)
val <- diag(exp(1:4)/(1+exp(1:4)))
})

test_that("Test domain for huber", {
value(a) <- 2
expr <- CVXR::huber(a)

value(a) <- 3
expr <- CVXR::huber(a, M = 2)

value(a) <- -1
expr <- CVXR::huber(a, M = 2)

value(x) <- c(3,4)
expr <- CVXR::huber(x)
val <- diag(c(2,2))

value(x) <- c(-1e-9,4)
expr <- CVXR::huber(x)
val <- diag(c(0,2))

value(A) <- cbind(c(1,2), c(3,4))
expr <- CVXR::huber(A, M = 3)
val <- diag(c(2,4,6,6))
})

test_that("Test domain for kl_div", {
b <- Variable()
value(a) <- 2
value(b) <- 4
expr <- kl_div(a, b)

value(a) <- 3
value(b) <- 0
expr <- kl_div(a, b)

value(a) <- -1
value(b) <- 2
expr <- kl_div(a, b)

y <- Variable(2)
value(x) <- c(3,4)
value(y) <- c(5,8)
expr <- kl_div(x, y)
val <- diag(log(c(3,4)) - log(c(5,8)))
val <- diag(c(1-3/5,1-4/8))

value(x) <- c(-1e-9,4)
value(y) <- c(1,2)
expr <- kl_div(x, y)

value(A) <- cbind(c(1,2), c(3,4))
value(B) <- cbind(c(5,1), c(3.5,2.3))
expr <- kl_div(A, B)
div <- as.numeric(value(A) / value(B))
val <- diag(log(div))
val <- diag(1-div)
})

test_that("Test domain for max_elemwise", {
b <- Variable()
value(a) <- 2
value(b) <- 4
expr <- max_elemwise(a, b)

value(a) <- 3
value(b) <- 0
expr <- max_elemwise(a, b)

value(a) <- -1
value(b) <- 2
expr <- max_elemwise(a, b)

y <- Variable(2)
value(x) <- c(3,4)
value(y) <- c(5,-5)
expr <- max_elemwise(x, y)
val <- diag(c(0,1))
val <- diag(c(1,0))

value(x) <- c(-1e-9,4)
value(y) <- c(1,4)
expr <- max_elemwise(x, y)
val <- diag(c(0,1))
val <- diag(c(1,0))

value(A) <- cbind(c(1,2), c(3,4))
value(B) <- cbind(c(5,1), c(3,2.3))
expr <- max_elemwise(A, B)
div <- as.numeric(value(A) / value(B))
val <- diag(c(0,1,1,1))
val <- diag(c(1,0,0,0))
})

test_that("Test domain for min_elemwise", {
b <- Variable()
value(a) <- 2
value(b) <- 4
expr <- min_elemwise(a, b)

value(a) <- 3
value(b) <- 0
expr <- min_elemwise(a, b)

value(a) <- -1
value(b) <- 2
expr <- min_elemwise(a, b)

y <- Variable(2)
value(x) <- c(3,4)
value(y) <- c(5,-5)
expr <- min_elemwise(x, y)
val <- diag(c(1,0))
val <- diag(c(0,1))

value(x) <- c(-1e-9,4)
value(y) <- c(1,4)
expr <- min_elemwise(x, y)
val <- diag(c(1,1))
val <- diag(c(0,0))

value(A) <- cbind(c(1,2), c(3,4))
value(B) <- cbind(c(5,1), c(3,2.3))
expr <- min_elemwise(A, B)
div <- as.numeric(value(A) / value(B))
val <- diag(c(1,0,1,0))
val <- diag(c(0,1,0,1))
})

test_that("Test domain for power", {
value(a) <- 2
expr <- sqrt(a)

value(a) <- 3
expr <- sqrt(a)

value(a) <- -1
expr <- sqrt(a)

value(x) <- c(3,4)
expr <- x^3

value(x) <- c(-1e-9,4)
expr <- x^3

value(A) <- cbind(c(1,-2), c(3,4))
expr <- A^2
val <- diag(c(2,-4,6,8))

# Constant
expr <- a^0

expr <- x^0
expect_equivalent(as.matrix(grad(expr)[[as.character(x@id)]]), matrix(0, nrow = 2, ncol = 2))
})

# test_that("Test grad for partial minimization/maximization problems", {
# for(obj in list(Minimize(a^-1), Maximize(Entr(a)))) {
#   prob <- Problem(obj, list(x + a >= c(5,8)))
#
#   # Optimize over nothing
#   expr <- partial_optimize(prob, dont_opt_vars = list(x, a))
#   value(a) <- NA
#   value(x) <- NA
#
#   # Outside domain
#   value(a) <- 1.0
#   value(x) <- c(5,5)
#
#   value(a) <- 1
#   value(x) <- c(10,10)
#
#   # Optimize over x
#   expr <- partial_optimize(prob, opt_vars = list(x))
#   value(a) <- 1
#
#   # Optimize over a
#   fix_prob <- Problem(obj, list(x + a >= c(5,8), x == 0))
#   result <- solve(fix_prob)
#   dual_val <- [email protected][1]@[email protected]
#   expr <- partial_optimize(prob, opt_vars = list(a))
#   value(x) <- c(0,0)
#
#   # Optimize over x and a
#   expr <- partial_optimize(prob, opt_vars = list(x, a))
# }
# })

test_that("Test grad for affine atoms", {
value(a) <- 2
expr <- -a

value(a) <- 2
expr <- 2*a

value(a) <- 2
expr <- a/2

value(x) <- c(3,4)
expr <- -x
val <- -diag(c(1,1))

value(A) <- cbind(c(1,2), c(3,4))
expr <- -A
val <- -diag(rep(1,4))

value(A) <- cbind(c(1,2), c(3,4))
expr <- A[1,2]
val <- matrix(0, nrow = 4, ncol = 1)
val[3] <- 1

z <- Variable(3)
value(x) <- c(1,2)
value(z) <- c(1,2,3)
expr <- vstack(x, z)
val <- matrix(0, nrow = 2, ncol = 5)
val[,1:2] <- diag(2)

val <- matrix(0, nrow = 3, ncol = 5)
val[,3:ncol(val)] <- diag(3)

# Cumulative sum
value(x) <- c(1,2)
expr <- cumsum_axis(x)
val <- matrix(1, nrow = 2, ncol = 2)
val[2,1] <- 0