在 Julia 中切割矩阵的更快方法

Faster way to slice a matrix in Julia

我想根据索引变量迭代不同的矩阵块。您可以将其视为如何在使用面板数据的模型上计算不同个体对对数似然的个体贡献。话虽如此,我希望它尽可能快。

我已经阅读了一些与之相关的问题。但是 none 他们直接回答了我的问题。例如,What is the recommended way to iterate a matrix over rows? shows ways to run over the WHOLE bunch of rows not iterating over blocks of rows. Additionally, Julia: loop over rows of matrix (or not) 也是关于如何再次遍历每一行而不是遍历它们的块。

所以这是我的问题。假设您有 X,这是一个 2x9 矩阵和一个索引样本中个体的 id 变量。我想迭代它们以尽快构建我的对数似然贡献。我在这里只是通过使用布尔值对矩阵进行切片来做到这一点,但这似乎相对低效,因为我 对每个人 检查整个向量以查看它们是否匹配。

id = [1 1 2 2 2 2 3 3 3 ]
x1 = [1 2 3 4 5 6 7 8 9]
x2 = [10 11 12 13 14 15 16 17 18]
X = transpose([x1 ; x2 ])
# 9×2 transpose(::Matrix{Int64}) with eltype Int64:
# 1  10
# 2  11
# 3  12
# 4  13
# 5  14
# 6  15
# 7  16
# 8  17
# 9  18

#Iteraring over the block of observations. 
for n in unique(id)
  select_row = vec(isone.(n.==transpose(id)))
  X_n = X[ select_row , : ]
  print(X_n)
end
#[1 10; 2 11]
#[3 12; 4 13; 5 14; 6 15]
#[7 16; 8 17; 9 18]

在其他语言中(见下文),我创建了一个矩阵,其中包含事先(在循环之外)每个人的位置,然后根据它索引我的回归量。不过,我不知道如何在 Julia 中执行此操作(或者是否有相应的函数),而且我也不知道这是否是使用这种语言的最快方法。


Mata 示例

例如,为了便于说明,在 Mata(Stata 的矩阵语言)中我会执行以下操作:

mata:
id = 1       \ 3\ 3 
x1 = 1       \ 8\ 9
x2 = 10       \ 17\ 18
X  = x1 , x2
/*this indicate where is each invididual*/
paninfo = panelsetup(id,1) 
paninfo
//       1   2
//    +---------+
//  1 |  1   2  |
//  2 |  3   6  |
//  3 |  7   9  |
//    +---------+
/*this gave me the total number of individuals*/
npanels = panelstats(paninfo)[1]  
npanels 
// 3
for(n=1; n <= npanels; ++n) {
    /* this selects the block of variables for individual n*/
    x_n = panelsubmatrix(X, n, paninfo)  
    x_n 
}
//        1    2
//    +-----------+
//  1 |   1   10  |
//  2 |   2   11  |
//    +-----------+
//        1    2
//    +-----------+
//  1 |   3   12  |
//  2 |   4   13  |
//  3 |   5   14  |
//  4 |   6   15  |
//    +-----------+
//        1    2
//    +-----------+
//  1 |   7   16  |
//  2 |   8   17  |
//  3 |   9   18  |
//    +-----------+
end 

首先,我建议您使用向量而不是矩阵:

julia> id = [1, 1, 2, 2, 2, 2, 3, 3, 3]
9-element Vector{Int64}:
 1
 1
 2
 2
 2
 2
 3
 3
 3

julia> x1 = [1, 2, 3, 4, 5, 6, 7, 8, 9]
9-element Vector{Int64}:
 1
 2
 3
 4
 5
 6
 7
 8
 9

julia> x2 = [10, 11, 12, 13, 14, 15, 16, 17, 18]
9-element Vector{Int64}:
 10
 11
 12
 13
 14
 15
 16
 17
 18

julia> X = [x1 x2]
9×2 Matrix{Int64}:
 1  10
 2  11
 3  12
 4  13
 5  14
 6  15
 7  16
 8  17
 9  18

现在,在您的上下文中对 id 进行 X 分组的相对简单的方法是使用 SplitApplyCombine.jl:

julia> res = map(x -> view(X, x, :), groupfind(id))
3-element Dictionaries.Dictionary{Int64, SubArray}
 1 │ [1 10; 2 11]
 2 │ [3 12; 4 13; 5 14; 6 15]
 3 │ [7 16; 8 17; 9 18]

julia> res[1]
2×2 view(::Matrix{Int64}, [1, 2], :) with eltype Int64:
 1  10
 2  11

julia> res[2]
4×2 view(::Matrix{Int64}, [3, 4, 5, 6], :) with eltype Int64:
 3  12
 4  13
 5  14
 6  15

julia> res[3]
3×2 view(::Matrix{Int64}, [7, 8, 9], :) with eltype Int64:
 7  16
 8  17
 9  18

我使用视图来提高操作的性能。