在Julia中创建一个稀疏对称随机矩阵

Question

有没有一种简单的方法在Julia中创建一个稀疏对称随机矩阵？Julia有命令“sprand(m，n，d)”，它“创建一个密度为d的m×n随机矩阵（稀疏），其中iid非零元素在半开区间[0,1)上均匀分布。”但据我所知，这并不一定返回对称矩阵。我正在寻找与MATLAB的“R = sprandsym(n，density)”相当的命令，它自动创建一个稀疏对称随机矩阵。如果还没有实现这样的命令，那么将“sprand(m，n，d)”返回的矩阵转换为对称矩阵的解决方法是什么？谢谢！

- miga89

Answer 1

You could Symmetric(sprand(10,10,0.4))

- Michael K. Borregaard

---

谢谢！生成的矩阵仍然是稀疏的吗？还是我需要执行 ´´´sparse(Symmetric(randn(m,n)))´´´ ？ - miga89

它仍然会很稀疏。对称矩阵类型包装其他矩阵并忽略它们的下三角部分。一个小注意点是，在底层的稀疏矩阵中存储无法访问的条目所使用的浪费内存。 - Dan Getz

{btsdaf} - Gnimuc

我提出的方法有一个优点，那就是它是通用的，并且将受益于未来对对称稀疏矩阵支持的改进。 - Michael K. Borregaard

Answer 2

为了避免Michael Borregaard回答评论中提到的额外内存限制，以下函数接受一个稀疏矩阵并删除下三角部分的条目。如果SparseMatrixCSC格式不熟悉，它还可以作为演示如何操作该格式的好方法：

function droplower(A::SparseMatrixCSC)
    m,n = size(A)
    rows = rowvals(A)
    vals = nonzeros(A)
    V = Vector{eltype(A)}()
    I = Vector{Int}()
    J = Vector{Int}()
    for i=1:n
        for j in nzrange(A,i)
            rows[j]>i && break
            push!(I,rows[j])
            push!(J,i)
            push!(V,vals[j])
        end
    end
    return sparse(I,J,V,m,n)
end

使用示例：

julia> a = [0.5 1.0 0.0 ; 2.0 0.0 0.0 ; 0.0 0.0 0.0]
3×3 Array{Float64,2}:
 0.5  1.0  0.0
 2.0  0.0  0.0
 0.0  0.0  0.0

julia> b = sparse(a)
3×3 SparseMatrixCSC{Float64,Int64} with 3 stored entries:
  [1, 1]  =  0.5
  [2, 1]  =  2.0
  [1, 2]  =  1.0

julia> c = droplower(b)
3×3 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 1]  =  0.5
  [1, 2]  =  1.0

julia> full(Symmetric(c))    # note this is symmetric although c isn't
3×3 Array{Float64,2}:
 0.5  1.0  0.0
 1.0  0.0  0.0
 0.0  0.0  0.0

- Dan Getz

---

is c symmetric? - Liso

@Liso c 不是对称的，但 Symmetric(c) 是对称的，并且它不会浪费存储空间在矩阵的下三角部分（正如我在第一段中所表达的）。 - Dan Getz

有triu（对应于您的droplower），以及tril来获取上三角或下三角。同样，还有triu!和tril!可以原地执行操作。 - fredrikekre

@fredrikekre 真不错！谢谢 - 有这样一个函数是很合理的。 - Dan Getz

Answer 3

对于SparseMatrixCSC上的操作，通常需要定制以获得最大效率。因此，为了从一个稀疏矩阵A转换为具有相同上部的对称稀疏矩阵，这里提供一个定制版本（它可能有点晦涩，但是有效）：

function symmetrize(A::SparseMatrixCSC)
    m,n = size(A)
    m == n || error("argument expected to be square matrix")
    rows = rowvals(A) ; vals = nonzeros(A)
    a = zeros(Int,n) ; b = zeros(Int,n) ; c = 0
    for i=1:n
        for j in nzrange(A, i)
            if rows[j]>=i
                if rows[j]==i a[i] += 1 ; c += 1 ; end
                break
            end
            a[i] += 1 ; b[rows[j]] += 1 ; c += 2
        end
    end
    c == 0 && return SparseMatrixCSC(n, n, ones(n+1), nrows, nvals)
    ncolptr = Vector{Int}(n+1)
    nrows = Vector{Int}(c) ; nvals = Vector{eltype(A)}(c)
    idx = 1
    for i=1:n
        ncolptr[i] = idx
        if a[i]==0 a[i] = idx ; idx += b[i] ; continue ; end
        for j in (0:a[i]-1)+first(nzrange(A, i))
            nvals[idx] = vals[j] ; nrows[idx] = rows[j] ; idx += 1
            rows[j] >= i && break
            nvals[a[rows[j]]] = vals[j] ; nrows[a[rows[j]]] = i
            a[rows[j]] += 1
        end
        a[i] = idx ; idx += b[i]
    end
    ncolptr[n+1] = idx
    return SparseMatrixCSC(n, n, ncolptr, nrows, nvals)
end

一个样例运行：

julia> f = sprand(5,5,0.2)
5×5 SparseMatrixCSC{Float64,Int64} with 5 stored entries:
  [1, 1]  =  0.981579
  [3, 1]  =  0.330961
  [5, 1]  =  0.527683
  [4, 5]  =  0.196898
  [5, 5]  =  0.579006

julia> full(f)
5×5 Array{Float64,2}:
 0.981579  0.0  0.0  0.0  0.0     
 0.0       0.0  0.0  0.0  0.0     
 0.330961  0.0  0.0  0.0  0.0     
 0.0       0.0  0.0  0.0  0.196898
 0.527683  0.0  0.0  0.0  0.579006

julia> full(symmetrize(f))
5×5 Array{Float64,2}:
 0.981579  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       0.196898
 0.0       0.0  0.0  0.196898  0.579006

这个版本应该比其他版本更快，但仍需要进行基准测试（并在for循环中添加一些@inbounds）。

- Dan Getz

---

你可以在Gnimuc链接的问题上对那段代码进行注释 - 它看起来像是他们正在考虑实现的那种东西。 - Michael K. Borregaard