使用逻辑索引进行稀疏矩阵分配的结果不佳

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

所以我得到了和提问者一样的东西。但是当我做事时 "my way"

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 一起工作。

一些额外的想法（编辑）

当我进一步思考这个问题时，我突然想到 sparse() 没有 BitArray 方法的一个原因是，为已经非常高的数据实现稀疏存储并不是很有用紧凑型。从上面考虑 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 方法也许不是疏忽。它代表显着 space 节省的情况可能不像人们乍看之下想象的那么普遍。

使用逻辑索引进行稀疏矩阵分配的结果不佳

Bad results in sparse matrix assignment with logical indexing

julia