干净的方式来计算雅可比阵的总和

Question

我在R做了一些优化，并与我需要编写一个函数，返回一个雅可比。这是一个非常简单的雅可比 - 只是零和一个 - 但我想快速，干净地填充它。我目前的代码工作，但很sl。。我有一个四维概率数组。通过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？

这是我的马虎溶液（我检查针对雅可比功能正确性从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作为输入 - 即，我的概率阵列的扁平表示 - - 因为优化函数需要向量参数。

Answer 1

花一些时间去了解你正在尝试做的，但这里是一个命题，以取代constraint_function_jacobian：1

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 
} 

除非我错了，雅可比建设填充对角线，因为它是不是我们必须将它分割上array_dim[4]方阵，以填补他们的对角线用方阵1.

我没有摆脱prob_vector转型的一个数组，然后获取其dim因为这将是一样的array_dim，跳过这一步是没有的这是一个巨大的改进，但它简化了IMO的代码。

结果根据测试都ok：

> 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