performance - 如何有效地构造一个大型 SparseArray？包这个？

Question

问题
Julia 有没有有效的方式 从给定的条目列表 (u,v,w) 构建一个巨大的稀疏矩阵 ，其中一些可以具有 相同地点 (u,v)，在这种情况下，它们的权重 w 必须是 求和 .因此u,v,w是输入向量，我希望创建一个具有值 w[i] 的稀疏矩阵在位置 u[i],v[i] .例如，Mathematica 代码

n=10^6;   m=500*n;   
u=RandomInteger[{1,n},m];   
v=RandomInteger[{1,n},m];   
w=RandomInteger[{-9, 9}, m];   AbsoluteTiming[
SetSystemOptions["SparseArrayOptions"->{"TreatRepeatedEntries"->1}]; 
a= SparseArray[{u,v}\[Transpose] -> w, {n,n}]; ]

需要 135 秒和 60GB 的 RAM。等效的 Python 代码

import scipy.sparse as sp   
import numpy as np
import time
def ti(): return time.perf_counter() 
n=10**6; m=500*n;
u=np.random.randint(0,n,size=m);            
v=np.random.randint(0,n,size=m);            
w=np.random.randint(-9,10,size=m); t0=ti(); 
a=sp.csr_matrix((w,(u,v)),shape=(n,n),dtype=int); t1=ti(); print(t1-t0)

需要 36 秒和 20GB，但不支持 (2)。等效的 Julia 代码

using SparseArrays;
m=n=10^6; r=500*n;   
u=rand(1:m,r);   
v=rand(1:n,r);   
w=rand(-9:9,r);
function f() a=sparse(u,v,w,m,n,+); end;   
@time f()

需要 155 秒和 26GB(@time 宏错误地报告它只使用了 15GB)。
有没有办法到 使这种结构更高效 ?
有没有 Julia 包 哪个更擅长这个？
背景
我为线性代数和同调代数计算创建了一个包。我在 Mathematica 中做到了，因为 SparseArray实现 (0) 有

非常高效(快速且低 RAM 使用率)，

支持精确分数作为矩阵条目。

支持多项式作为矩阵项。

然而，它不是

并行化

开源

启用 GPU 或集群。

出于各种长期原因，我正在考虑迁移到不同的编程语言。
我试过 Python SciPy.sparse ，但不满足(2)。我不确定 C++ 库是否支持 (2)。作为测试，我尝试对大小为 10^9 的浮点数组进行排序，并比较 Mathematica、Python numpy 的性能。和 Julia ThreadsX .后者给我留下了难以置信的印象(比 numpy 快 3 倍，比 Mathematica 快 10 倍)。我正在考虑迁移到 Julia，但我希望首先确保我的库会比在 Mathematica 中表现得更好。
附言我怎样才能实现我的 f()以上在调用时不打印/输出任何内容？
P.P.S.我也问过这个问题 here .

Answer 1

在这种特殊情况下，您似乎可以更简单地使用(例如)

julia> @time a = sprand(Float64, 10^6, 10^6, 500/10^6)
 17.880987 seconds (7 allocations: 7.459 GiB, 7.71% gc time)
1000000×1000000 SparseMatrixCSC{Float64, Int64} with 499998416 stored entries:
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛⠛

更一般地说，您可能想查看底层方法 SparseArrays.sparse! ，这将允许更有效的就地施工。

help?> SparseArrays.sparse!
  sparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv},
          m::Integer, n::Integer, combine, klasttouch::Vector{Ti},
          csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, csrnzval::Vector{Tv},
          [csccolptr::Vector{Ti}], [cscrowval::Vector{Ti}, cscnzval::Vector{Tv}] ) where {Tv,Ti<:Integer}

  Parent of and expert driver for sparse; see sparse for basic usage. This method allows
  the user to provide preallocated storage for sparse's intermediate objects and result
  as described below. This capability enables more efficient successive construction of
  SparseMatrixCSCs from coordinate representations, and also enables extraction of an
  unsorted-column representation of the result's transpose at no additional cost.

  This method consists of three major steps: (1) Counting-sort the provided coordinate
  representation into an unsorted-row CSR form including repeated entries. (2) Sweep
  through the CSR form, simultaneously calculating the desired CSC form's column-pointer
  array, detecting repeated entries, and repacking the CSR form with repeated entries
  combined; this stage yields an unsorted-row CSR form with no repeated entries. (3)
  Counting-sort the preceding CSR form into a fully-sorted CSC form with no repeated
  entries.

  Input arrays csrrowptr, csrcolval, and csrnzval constitute storage for the
  intermediate CSR forms and require length(csrrowptr) >= m + 1, length(csrcolval) >=
  length(I), and length(csrnzval >= length(I)). Input array klasttouch, workspace for
  the second stage, requires length(klasttouch) >= n. Optional input arrays csccolptr,
  cscrowval, and cscnzval constitute storage for the returned CSC form S. csccolptr
  requires length(csccolptr) >= n + 1. If necessary, cscrowval and cscnzval are
  automatically resized to satisfy length(cscrowval) >= nnz(S) and length(cscnzval) >=
  nnz(S); hence, if nnz(S) is unknown at the outset, passing in empty vectors of the
  appropriate type (Vector{Ti}() and Vector{Tv}() respectively) suffices, or calling the
  sparse! method neglecting cscrowval and cscnzval.

  On return, csrrowptr, csrcolval, and csrnzval contain an unsorted-column
  representation of the result's transpose.

  You may reuse the input arrays' storage (I, J, V) for the output arrays (csccolptr,
  cscrowval, cscnzval). For example, you may call sparse!(I, J, V, csrrowptr, csrcolval,
  csrnzval, I, J, V).

  For the sake of efficiency, this method performs no argument checking beyond 1 <= I[k]
  <= m and 1 <= J[k] <= n. Use with care. Testing with --check-bounds=yes is wise.

  This method runs in O(m, n, length(I)) time. The HALFPERM algorithm described in F.
  Gustavson, "Two fast algorithms for sparse matrices: multiplication and permuted
  transposition," ACM TOMS 4(3), 250-269 (1978) inspired this method's use of a pair of
  counting sorts.

如果您需要比这更快的东西，原则上您可以使用 SparseMatrixCSC 构建一个几乎没有开销的稀疏矩阵。直接构造函数(尽管您必须知道自己在做什么以及字段的含义)。作为现在 SparseArrays 标准库的一些早期贡献者(在 Base 的时间部分) noted

I consider using the SparseMatrixCSC constructor directly to be the I-know-what-I'm-doing interface. People use it to skip the overhead and checking of sparse, for example to create matrices with non-sorted rows, or explicit zeros, things like that. I'd be fine with making show do some validity checks on the assumptions it makes before trying to show the matrix, but having a lenient constructor as a direct way of saying "trust me, wrap this data in a SparseMatrixCSC" is a valuable thing.

如果你愿意，你当然可以编写并行化函数(使用 ThreadsX.jl 或 LoopVectorization.jl 的线程版本或其他任何东西)来直接并行构造底层数组，然后用 SparseMatrixCSC 包裹它们。
有关如何构建 SparseMatrixCSC 的所有详细信息您可以直接查看 https://github.com/JuliaLang/julia/blob/master/stdlib/SparseArrays/src/sparsematrix.jl
编辑:举一个说明性的例子(虽然见下面的警告)

using SparseArrays
m = n = 10^6
colptr = collect(1:500:n*500+1)
rowval = repeat(1:2000:m, n)
nzval = rand(Float64, n*500)
@time a = SparseMatrixCSC(m,n,colptr,rowval,nzval)

julia> @time a = SparseMatrixCSC(m,n,colptr,rowval,nzval)
  0.012364 seconds (1 allocation: 48 bytes)
1000000×1000000 SparseMatrixCSC

{Float64, Int64} with 500000000 stored entries:
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿

julia> a[1,1]
0.9356478932120766

julia> a[2,1]
0.0

julia> a[2001,1]
0.2121877970335444

julia> a[:,1]
1000000-element SparseVector{Float64, Int64} with 500 stored entries:
  [1      ]  =  0.935648
  [2001   ]  =  0.212188
  [4001   ]  =  0.429638
  [6001   ]  =  0.0190535
  [8001   ]  =  0.0878085
  [10001  ]  =  0.24005
  [12001  ]  =  0.785151
  [14001  ]  =  0.348142
             ⋮
  [982001 ]  =  0.637904
  [984001 ]  =  0.136397
  [986001 ]  =  0.483078
  [988001 ]  =  0.862434
  [990001 ]  =  0.703863
  [992001 ]  =  0.00990588
  [994001 ]  =  0.629455
  [996001 ]  =  0.123507
  [998001 ]  =  0.411807

两个主要警告:

在此示例中没有检查行是否已排序。如果您不知道这意味着什么或为什么这很重要，那么这种低级方法不适合您。

这不提供冲突检查。如果您已经从数据来源知道没有冲突，那很好。对于随机矩阵，它在统计上可能好也可能不好，这取决于稀疏程度。但同样，如果您不知道这意味着，请远离。