乾淨的方式來計算雅可比陣的總和

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