# Populate a symmetric sparse matrix

**URL:** <https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830>\
**Category:** General Usage\
**Tags:** question, sparse\
**Created:** [August 15, 2023, 2:58pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830 "2023-08-15T14:58:45Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 15, 2023, 2:58pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/1 "2023-08-15T14:58:45Z")

</div>

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:

```julia
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,

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 15, 2023, 3:35pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/2 "2023-08-15T15:35:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [August 15, 2023, 3:45pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/3 "2023-08-15T15:45:56Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 15, 2023, 4:39pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/4 "2023-08-15T16:39:24Z")

</div>

> [@Dan](#):
>
> ```julia
> for i in 1:order
> for j in i:order
> sparsity_pattern(i,j,order) || continue
> 
> ```

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.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [August 15, 2023, 4:51pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/5 "2023-08-15T16:51:49Z")

</div>

> [@stevengj](#):
>
> No, this is still O(n^2) for an n \times nn×nn \times n matrix

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.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 15, 2023, 5:30pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/6 "2023-08-15T17:30:08Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 15, 2023, 6:50pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/7 "2023-08-15T18:50:57Z")

</div>

> [@mleprovost](#):
>
> For a tuple of integers …

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\%).

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 15, 2023, 7:01pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/8 "2023-08-15T19:01:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 15, 2023, 7:03pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/9 "2023-08-15T19:03:20Z")

</div>

> [@mleprovost](#):
>
> Are there techniques to determine the pattern of nonzeros entries from this condition for a fixed m?

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

```julia
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).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 15, 2023, 7:07pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/10 "2023-08-15T19:07:09Z")

</div>

> [@mleprovost](#):
>
> Thank you for you answer. In my numerical experiments, the fill-in fraction for a given m is usually about 20%.

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
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
julia> @btime $Asparse * $Asparse;
  51.861 ms (6 allocations: 15.27 MiB)

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

```

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [August 15, 2023, 7:16pm UTC](https://discourse.julialang.org/t/populate-a-symmetric-sparse-matrix/102830/11 "2023-08-15T19:16:42Z")

</div>

Thanks a lot for your insights.