julia - 带有逻辑索引的稀疏矩阵分配的不良结果

Question

在 Matlab/Octave 中，我可以使用逻辑索引在矩阵 A 中满足特定要求的每个位置为矩阵 B 分配一个值。

octave:1> A = [.1;.2;.3;.4;.11;.13;.14;.01;.04;.09];
octave:2> C = A < .12
C =

   1
   0
   0
   0
   1
   0
   0
   1
   1
   1

octave:3> B = spalloc(10,1);
octave:4> B(C) = 1
B =

Compressed Column Sparse (rows = 10, cols = 1, nnz = 5 [50%])

  (1, 1) ->  1
  (5, 1) ->  1
  (8, 1) ->  1
  (9, 1) ->  1
  (10, 1) ->  1

但是，如果我在 Julia 中尝试基本相同的代码，结果是不正确的:

julia> A = [.1;.2;.3;.4;.11;.13;.14;.01;.04;.09];

julia> B = spzeros(10,1)
10x1 sparse matrix with 0 Float64 entries:

julia> C = A .< .12
10-element BitArray{1}:
  true
 false
 false
 false
  true
 false
 false
  true
  true
  true

julia> B[C] = 1
1

julia> B
10x1 sparse matrix with 5 Float64 entries:
    [0 ,  1]  =  1.0
    [0 ,  1]  =  1.0
    [1 ,  1]  =  1.0
    [1 ,  1]  =  1.0
    [1 ,  1]  =  1.0

我是否在某个地方的语法中犯了错误，是我误解了什么，还是这是一个错误？注意，如果我在 Julia 中使用完整矩阵，我会得到正确的结果，但由于我的应用程序中的矩阵非常稀疏(有限元模拟中的基本边界条件)，我更喜欢使用稀疏矩阵

Answer 1

看起来好像sparse BitArray有一些问题的。

julia> VERSION
v"0.3.5"
julia> A = [.1;.2;.3;.4;.11;.13;.14;.01;.04;.09]
julia> B = spzeros(10,1)
julia> C = A .< .12
julia> B[C] = 1
julia> B
10x1 sparse matrix with 5 Float64 entries:
        [0 ,  1]  =  1.0
        [0 ,  1]  =  1.0
        [1 ,  1]  =  1.0
        [1 ,  1]  =  1.0
        [1 ,  1]  =  1.0

所以我得到了和提问者一样的东西。然而，当我以“我的方式”做事时

julia> B = sparse(C)
ERROR: `sparse` has no method matching sparse(::BitArray{1})

julia> B = sparse(float(C))
10x1 sparse matrix with 5 Float64 entries:
        [1 ,  1]  =  1.0
        [5 ,  1]  =  1.0
        [8 ,  1]  =  1.0
        [9 ,  1]  =  1.0
        [10,  1]  =  1.0

因此，如果您转换 BitArray，这将起作用至 Float .我想这个解决方法会让你继续前进，但似乎 sparse应该与 BitArray 一起使用.

一些额外的想法(编辑)

当我进一步思考这一点时，我想到了没有 BitArray 的一个原因。 sparse() 的方法是为已经高度紧凑的类型实现稀疏存储并不是非常有用。考虑 B和 C从上面:

julia> sizeof(C)
8

julia> sizeof(B)
40

所以对于这些数据， sparse版本比原版大很多。乍一看，它实际上比这个简单(也许过于简单)的检查显示更糟糕。 sizeof(::BitArray{1})看起来是整个数组的大小，但是 sizeof(::SparseMatrixCSC{})显示存储的每个元素的大小。所以真正的大小差异是 8 字节和 200 字节。

当然，如果数据足够稀疏(略低于 1% true )，稀疏存储开始胜出，尽管它的开销很高。

julia> C = rand(10^6) .< 0.01

julia> B = sparse(float(C))

julia> sizeof(C)
125000

julia> sum(C)*sizeof(B)
394520

julia> C = rand(10^6) .< 0.001

julia> B = sparse(float(C))

julia> sizeof(C)
125000

julia> sum(C)*sizeof(B)
40280

因此， sparse() 或许不是疏忽。没有 BitArray方法。它代表显着节省空间的情况可能不像乍一看想象的那么常见。