# Sparse syrk() routine using Transducers.jl

**URL:** <https://discourse.julialang.org/t/sparse-syrk-routine-using-transducers-jl/86218>\
**Category:** Performance\
**Tags:** blas, sparse, transducers\
**Created:** [August 23, 2022, 7:09pm UTC](https://discourse.julialang.org/t/sparse-syrk-routine-using-transducers-jl/86218 "2022-08-23T19:09:07Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [August 23, 2022, 7:09pm UTC](https://discourse.julialang.org/t/sparse-syrk-routine-using-transducers-jl/86218/1 "2022-08-23T19:09:07Z")

</div>

Hi, I implemented a sparse [syrk](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.syrk) routine to see if I could get something faster than the naive implementation:

```julia

using SparseArrays
using LinearAlgebra
using Transducers
using Combinatorics

naive_spsyrk(s) = Symmetric(s * s')

spsyrk(s) =
    zip(findnz(s)...) |>
    PartitionBy(x -> x[2]) |>
    MapCat(a -> with_replacement_combinations(a, 2)) |>
    Map(x -> (x[1][1], x[2][1], x[1][3] * x[2][3])) |>
    foldxl(TeeRF(
               Map(x -> x[1])'(push!!),
               Map(x -> x[2])'(push!!),
               Map(x -> x[3])'(push!!),
           )) |>
    x -> sparse(x...) |>
    Symmetric

```

Unsurprisingly this attempt was unsuccessful, although I did get to within a less than an order of magnitude:

```julia
julia> a = sprand(1000, 1000, 1e-3);

julia> @benchmark naive_spsyrk($a)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 22.148 μs … 1.286 ms ┊ GC (min … max): 0.00% … 95.01%
 Time (median): 23.958 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 27.254 μs ± 45.658 μs ┊ GC (mean ± σ): 8.22% ± 4.83%

  ▂▆█▇▄▃▂▁ ▂
  █████████▅▆█▇▅▄▁▄▄▃▃▁▁▁▃▁▁▁▁▁▄▁▅▅▅▄▃▁▁▃▁▁▁▁▁▁▁▁▁▃▅▅▆▆▆▆▆▆▇▇ █
  22.1 μs Histogram: log(frequency) by time 61.6 μs <

 Memory estimate: 130.81 KiB, allocs estimate: 14.

julia> @benchmark spsyrk($a)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 166.617 μs … 4.131 ms ┊ GC (min … max): 0.00% … 93.70%
 Time (median): 173.007 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 214.630 μs ± 314.415 μs ┊ GC (mean ± σ): 15.76% ± 10.04%

  █▃▂ ▁
  ████▇▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ █
  167 μs Histogram: log(frequency) by time 2.92 ms <

 Memory estimate: 547.00 KiB, allocs estimate: 5763.

```

But maybe some folks here have ideas on speeding this implementation up. (My personal guess is the `PartitionBy` is the culprit here…)