Чистый способ вычисления якобиана суммирования массива

Я занимаюсь оптимизацией в R и в связи с этим мне нужно написать функцию, которая возвращает якобиан. Это очень простой якобиан — только нули и единицы — но я хотел бы заполнить его быстро и чисто. Мой текущий код работает, но очень неаккуратно.

У меня есть четырехмерный массив вероятностей. Проиндексируйте размеры с помощью i, j, k, l. Мое ограничение состоит в том, что для каждого i, j, k сумма вероятностей по индексу l должна равняться 1.

Я вычисляю свой вектор ограничений следующим образом:

get_prob_array_from_vector <- function(prob_vector, array_dim) {
    return(array(prob_vector, array_dim))
}

constraint_function <- function(prob_vector, array_dim) {
    prob_array <- get_prob_array_from_vector(prob_vector, array_dim)
    prob_array_sums <- apply(prob_array, MARGIN=c(1, 2, 3), FUN=sum)
    return(as.vector(prob_array_sums) - 1)  # Should equal zero
}

Мой вопрос: каков чистый и быстрый способ вычисления якобиана для as.vector(apply(array(my_input_vector, array_dim), MARGIN=c(1, 2, 3), FUN=sum)) -- т. е. моего constraint_function в приведенном выше коде -- по отношению к my_input_vector?

Вот мое небрежное решение (которое я проверяю на правильность с помощью функции jacobian из пакета numDeriv):

library(numDeriv)

array_dim <- c(5, 4, 3, 3)

get_prob_array_from_vector <- function(prob_vector, array_dim) {
    return(array(prob_vector, array_dim))
}

constraint_function <- function(prob_vector, array_dim) {
    prob_array <- get_prob_array_from_vector(prob_vector, array_dim)
    prob_array_sums <- apply(prob_array, MARGIN=c(1, 2, 3), FUN=sum)
    return(as.vector(prob_array_sums) - 1)
}

constraint_function_jacobian <- function(prob_vector, array_dim) {
    prob_array <- get_prob_array_from_vector(prob_vector, array_dim)
    jacobian <- matrix(0, Reduce("*", dim(prob_array)[1:3]), length(prob_vector))
    ## Must be a faster, clearner way of populating jacobian
    for(i in seq_along(prob_vector)) {
        dummy_vector <- rep(0, length(prob_vector))
        dummy_vector[i] <- 1
        dummy_array <- get_prob_array_from_vector(dummy_vector, array_dim)
        dummy_array_sums <- apply(dummy_array, MARGIN=c(1, 2, 3), FUN=sum)
        jacobian_row_idx <- which(dummy_array_sums != 0, arr.ind=FALSE)
        stopifnot(length(jacobian_row_idx) == 1)
        jacobian[jacobian_row_idx, i] <- 1
    }  # Is there a fast, readable one-liner that does the same as this for loop?
    stopifnot(sum(jacobian) == length(prob_vector))
    stopifnot(all(jacobian == 0 | jacobian == 1))
    return(jacobian)
}

## Example of a probability array satisfying my constraint
my_prob_array <- array(0, array_dim)
for(i in seq_len(array_dim[1])) {
    for(j in seq_len(array_dim[2])) {
        my_prob_array[i, j, , ] <- diag(array_dim[3])
    }
}
my_prob_array[1, 1, , ] <- 1 / array_dim[3]
my_prob_array[2, 1, , ] <- 0.25 * (1 / array_dim[3]) + 0.75 * diag(array_dim[3])

my_prob_vector <- as.vector(my_prob_array)  # Flattened representation of my_prob_array
should_be_zero_vector <- constraint_function(my_prob_vector, array_dim)
is.vector(should_be_zero_vector)
all(should_be_zero_vector == 0)  # Constraint is satistied

## Check constraint_function_jacobian for correctness using numDeriv
jacobian_analytical <- constraint_function_jacobian(my_prob_vector, array_dim)
jacobian_numerical <- jacobian(constraint_function, my_prob_vector, array_dim=array_dim)
max(abs(jacobian_analytical - jacobian_numerical))  # Very small

Мои функции принимают prob_vector в качестве входных данных, т. е. уплощенное представление моего массива вероятностей, потому что функциям оптимизации требуются векторные аргументы.


r arrays matrix-indexing derivative
person Adrian    schedule 04.03.2016    source источник

Ответы (1)


arrow_upward
7
arrow_downward

Потратьте некоторое время, чтобы понять, что вы пытались сделать, но вот предложение заменить ваш constraint_function_jacobian:

enhanced <- function(prob_vector,array_dim) {
  firstdim <- Reduce("*", array_dim[1:3])
  seconddim <- length(prob_vector)
  jacobian <- matrix(0, firstdim, seconddim)
  idxs <- split(1:seconddim,cut(1:seconddim,array_dim[4],labels=F))
  for( i in seq_along(idxs)) {
    diag(jacobian[, idxs[[i]] ]) <- 1
  }
  stopifnot(sum(jacobian) == length(prob_vector))
  stopifnot(all(jacobian == 0 | jacobian == 1))
  jacobian
}

Если я не ошибаюсь, якобианская конструкция заполняет диагонали 1, поскольку это не квадратная матрица, мы должны разделить ее на array_dim[4] квадратную матрицу, чтобы заполнить их диагонали 1.

Я избавился от преобразования prob_vector в массив, чтобы затем получить его dim, поскольку он будет таким же, как array_dim, пропуск этого шага не является огромным улучшением, но упрощает код IMO.

Результаты в порядке согласно тесту:

> identical(constraint_function_jacobian(my_prob_vector,array_dim),enhanced(my_prob_vector,array_dim))
[1] TRUE

Согласно бенчмарку, это дает большое ускорение:

> microbenchmark(constraint_function_jacobian(my_prob_vector,array_dim),enhanced(my_prob_vector,array_dim),times=100)

Unit: microseconds
                                                    expr       min        lq      mean     median         uq       max neval cld
 constraint_function_jacobian(my_prob_vector, array_dim) 16946.979 18466.491 20150.304 19066.7410 19671.4100 28148.035   100   b
                     enhanced(my_prob_vector, array_dim)   678.222   737.948   799.005   796.3905   834.5925  1141.773   100  a
person Tensibai    schedule 04.03.2016
comment
Спасибо, это серьезное улучшение! - person Adrian; 04.03.2016