julia - 在 Julia 中创建 UpperTriangularMatrix 的子矩阵

Question

给定一个列和行的列表，我想生成一个 cholesky 分解的子矩阵。示例:

julia>  A = rand(10,10)
julia>  R = chol(A'*A)
julia>  ind = [1,3,6,8,9]
julia>  R[ind,ind]

但是，这会导致错误:

ERROR: BoundsError: attempt to access 5x5
UpperTriangular{Float64,Array{Float64,2}}:
 1.28259  0.0           0.0  0.0  0.0
 0.0      6.51646e-314  0.0  0.0  0.0
 0.0      0.0           0.0  0.0  0.0
 0.0      0.0           0.0  0.0  0.0
 0.0      0.0           0.0  0.0  0.0
at index [2,1]
in _unsafe_getindex at multidimensional.jl:197
in getindex at abstractarray.jl:483

我知道这适用于典型的矩阵，但是 UpperTriangular 类型显然需要一些不同的东西......我找不到这方面的文档。

Answer 1

在 0.4 中，三角矩阵似乎没有更新以利用后备非标量索引(这是 0.3 中的一个缺失方法错误)。

目前最简单的解决方法是在索引之前转换为完整数组:

julia> full(R)[ind,ind]
5x5 Array{Float64,2}:
 2.2261  1.28096  1.69087    1.26135    1.50703
 0.0     1.03681  0.115735   0.559855   0.70766
 0.0     0.0      0.702936  -0.111155  -0.61263
 0.0     0.0      0.0        0.661491   0.33661
 0.0     0.0      0.0        0.0        0.159691

或者通过使用 SubArray，创建原始数据的 View (因此修改将传播):

julia> sub(R, ind, ind)
5x5 SubArray{Float64,2,UpperTriangular{Float64,Array{Float64,2}},Tuple{Array{Int64,1},Array{Int64,1}},0}:
 2.2261  1.28096  1.69087    1.26135    1.50703
 0.0     1.03681  0.115735   0.559855   0.70766
 0.0     0.0      0.702936  -0.111155  -0.61263
 0.0     0.0      0.0        0.661491   0.33661
 0.0     0.0      0.0        0.0        0.159691