Populate a symmetric sparse matrix

General Usage
,
1

Hello,

I would like to populate a sparse matrix A with a function f:\mathbb{N}\times \mathbb{N} \to \mathbb{R}, (i_1, i_2) \mapsto f(i_1, i_2) which is symmetric, i.e. f(i_1, i_2) = f(i_2, i_1). We want to fill A such that A[i_1, i_2] = f(i_1, i_2). The sparsity pattern of A is known and symmetric. I would like to reduce functions evaluations by using the symmetry of f when filling A.

If A is stored in the CSC format, how do you determine which indices correspond to a pair (i_1, i_2) and its symmetric (i_2, i_1) from A.colptr and A.rowval?

The sparsity pattern of A is given by a function \operatorname{sparsity\_pattern}(i_1, i_2). I would like to compute the sparsity pattern store it in a matrix pattern and transfer it when initializing A. Is there a better way to proceed?

Some MWE:

using LinearAlgebra
using SparseArrays 

f(i1, i2) = 2.0*(i1+i2)

# sparsity_pattern depends on a third entry i3 but is symmetric in i1 and i2
sparsity_pattern(i1, i2, i3) = (mod(i1 + i2 + i3, 2) == 0 && i1 + i2 >= i3 && i1 + i3 >= i2 && i2 + i3 >= i1)

order = 5

pattern = zeros(Int64, order, order)

for i1=1:order
    # Use symmetry
    for i2=i1:order
        if sparsity_pattern(i1, i2, order) == true 
            pattern[i1, i2] = 1
            pattern[i2, i1] = 1
        end
    end
end

pattern = sparse(pattern)
rows, cols, _ = findnz(pattern)

# Initialize A by transferring sparsity pattern of pattern
# Probably a better way using similar(pattern, Float64), not sure of the precise syntax
# A = SparseMatrixCSC(order, order, pattern.colptr, pattern.rowval, 
#                    Array{Float64}(undef,length(pattern.rowval)))

# Cleaner way
A = sparse(rows, cols, Array{Float64}(undef,length(rows)))

# Populate A without using the symmetry of f(i1, i2)

for (i1, i2) in zip(rows, cols)
    A[i1,i2] = f(i1, i2)
end

Thank you for your help,

2

Forget about the symmetry for now, you have deeper problems.

First, a 5x5 matrix is not big enough for a sparse-matrix data structure to be worthwhile. And in your example code your matrix is about 30% nonzero, which is not sparse enough to be worthwhile even if the matrix is large. How big and how sparse are your actual matrices?

Second, even if your actual matrix is much bigger and much sparser, you don’t want to allocate a dense matrix, loop over all the entries, and then convert to sparse — the goal in sparse-matrix computations is to have the work and storage be proportional only to the number of nonzero entries. So your fundamental approach of a sparsity_pattern function is wrongheaded.

3 Likes
3

Given the issues stevengj raised. An example of how to get around some of the problems:

R = Int[]; C = Int[]; V = Float64[]
for i in 1:order
    for j in i:order
        sparsity_pattern(i,j,order) || continue
        v = f(i,j)
        push!(R, i); push!(C, j); push!(V, v)
        push!(R, j); push!(C, i); push!(V, v)
    end
end
A = sparse(R, C, V, order, order)

In real usage, this should be in a function (and not global scope). A bit more optimizition can be done with reducing the number of push! calls, and if the counts of the nonzero elements are known in advance, pre-allocation could save time too.

1 Like
4

No, this is still O(n^2)` for an n \times n matrix. Querying every element to check whether it is nonzero does not scale.

You really have to figure out the pattern of nonzeros and just loop over those, so that your cost (in both time and space!) is proportional to the number of nonzeros only.

5

Right. O(n^2) for sparsity_pattern but ~nnz/2 for f function which could be the heavy part of the calculation. This version has nnz memory requirement.
In any case, the OP’s exact situation matters.

6

Thank you for sharing your insights.

For a tuple of integers (k,l,m), the entry (k,l) where we fix m is non-zero if k+l+m is even (\operatorname{mod}(k+l+m,2) = 0)) and the sum of any two integers is not smaller that the third one (k + l \geq m and k + m \geq l and l + m \geq k).

Are there techniques to determine the pattern of nonzeros entries from this condition for a fixed m?

7

This (corresponding to your sparsity_pattern formula above), is not sparse enough to benefit from a generic sparse-matrix data structure ala sparse, I think. You’re almost certainly better off just using a dense matrix (unless you want to develop a specialized “structured sparse” format for your application).

Dense matrices benefit from extremely regular memory-access patterns that can be (and are) highly optimized. Generic sparse matrices (e.g. the CSR format used by Julia) use a more complicated data structure, and are usually only more efficient if only a tiny fraction of the matrix is nonzero (usually \ll 10\%).

8

Thank you for you answer. In my numerical experiments, the fill-in fraction for a given m is usually about 20\%.

In my problem, the sparsity-pattern is given by element-wise product of these matrices for multiple values of the parameter m. For two values m and m^\prime, the sparsity pattern of the associated matrices is usually distinct, such that more sparsity is built after several multiplications of these matrices.

9

But yes, there are techniques, known as “algebra” :wink: . If you look at your matrix and think about the formulas you’ll see it’s quite a simple pattern:

julia> pattern = sparse(pattern)
5Ă—5 SparseMatrixCSC{Int64, Int64} with 10 stored entries:
 â‹…  â‹…  â‹…  1  â‹…
 â‹…  â‹…  1  â‹…  1
 â‹…  1  â‹…  1  â‹…
 1  â‹…  1  â‹…  1
 â‹…  1  â‹…  1  â‹…

If you think about it, you should be able to derive a simple formula for the first nonzero row in each column. And subsequently it is just every other row.

But if you only have 20% sparsity it might not be worth the bother. Just use something like @Dan’s solution (if you want to use sparse matrices at all).

1 Like
10

Depending on the size of your matrix, this might not be enough for sparsity to be worthwhile (again, with a generic sparse data structure, not something specialized to your sparsity pattern). e.g. for a 1000 \times 1000 matrix with 20% sparsity, on my machine the dense matrix–vector multiply is faster:

julia> m = 1000; Asparse = sprand(m,m,0.2); Adense = Matrix(Asparse);

julia> x = rand(m); y = similar(x);

julia> @btime mul!($y, $Asparse, $x);
  124.586 ÎĽs (0 allocations: 0 bytes)

julia> @btime mul!($y, $Adense, $x);
  80.081 ÎĽs (0 allocations: 0 bytes)

and the improvement is even bigger for matrix–matrix multiplies:

julia> @btime $Asparse * $Asparse;
  51.861 ms (6 allocations: 15.27 MiB)

julia> @btime $Adense * $Adense;
  13.697 ms (2 allocations: 7.63 MiB)
1 Like
11

Thanks a lot for your insights.