# Sum(log.(p \* C)) is 2 to 4 times slower than in NumPy

**URL:** <https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960>\
**Category:** Performance\
**Tags:** python, matrices\
**Created:** [January 26, 2021, 8:38am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960 "2021-01-26T08:38:50Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [January 26, 2021, 8:38am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/1 "2021-01-26T08:38:51Z")

</div>

I’m coming from Python and trying to port my Python code to Julia to gain more performance.

I want to compute a lot of log-likelihoods in a loop, so I need a fast function that computes the log-likelihood. The computation is essentially `sum(log.(p * C))`, where `C` is a constant matrix and `p` is a row vector.

Here’s my Python code:

```julia-auto
import numpy as np

def log_likelihood(p, C):
	return np.log(p @ C).sum()

Nx = 5500
Np = 200

p = np.linspace(0, 1000, Np)
C = np.random.rand(Np, Nx)

print(f"Python with Nx={Nx}, Np={Np}")
%timeit log_likelihood(p, C)

```

And Julia code:

```julia-auto
using BenchmarkTools
using LinearAlgebra

log_likelihood_types(p::Array{Float64, 2}, C::Array{Float64, 2})::Float64 = sum(log.(p * C))
log_likelihood_no_types(p, C) = sum(log.(p * C))

Nx = 5500
Np = 200

p = range(0, 1000, length=Np)' |> collect
C = rand(Np, Nx)

# Compile
log_likelihood_types(p, C)
log_likelihood_no_types(p, C)

println("Julia with Nx=$Nx, Np=$Np; BLAS threads: $(BLAS.get_num_threads())")

@btime log_likelihood_types($p, $C)
@btime log_likelihood_no_types($p, $C)

```

## Test results

```julia-auto
forcebru ~/t/julia_numpy_bench> ipython row_mul.ipy
Python with Nx=5500, Np=200
349 µs ± 1.45 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
forcebru ~/t/julia_numpy_bench> julia-1.6 row_mul.jl
Julia with Nx=5500, Np=200; BLAS threads: 4
  1.240 ms (4 allocations: 86.16 KiB)
  1.238 ms (4 allocations: 86.16 KiB)

```

So, 349 µs with NumPy and 1238 µs with Julia, so Julia is 3.5 times (!) slower. Also see this benchmark ([NumPy vs Julia · GitHub](https://gist.github.com/ForceBru/2da87f80725d62eb8f62454d53e00118)) which shows similar results.

For `Nx == 500` I’m getting:

- Julia: 101.575 μs
- NumPy: 36.9 µs ± 152 ns

Here Julia is 2.7 times slower!

## Versions

- julia version 1.6.0-beta1
- `np. __version__ == '1.19.5'`
- Python 3.7.1

## Questions

Am I doing something wrong? Am I maybe measuring something incorrectly? Why the 3x speed difference? Is there a way to “fix” it?

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [January 26, 2021, 8:58am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/2 "2021-01-26T08:58:53Z")

</div>

Better performance can be achieved by using

```julia
p = collect(range(0, 1000, length=Np))'

```

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [January 26, 2021, 8:59am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/3 "2021-01-26T08:59:10Z")

</div>

What makes a big difference is to keep your row vector a row vector type, not a Matrix:

```julia
julia> @btime sum(log.($p * $C));
  717.158 μs (4 allocations: 86.16 KiB)

julia> pa = collect(range(0, 1000, length=Np))' # type Adjoint, not 1xN Matrix
1×200 adjoint(::Vector{Float64}) with eltype Float64:
 0.0 5.02513 10.0503 15.0754 20.1005 25.1256 … 979.899 984.925 989.95 994.975 1000.0

julia> @btime sum(log.($pa * $C));
  239.876 μs (4 allocations: 86.16 KiB)

julia> @btime sum(log, $pa * $C);
  234.405 μs (2 allocations: 43.08 KiB)

julia> BLAS.vendor()
:openblas64

```

Same computer:

```julia
Python with Nx=5500, Np=200
236 µs ± 8.65 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

```

Using `sum(log, x)` saves allocating another array for the result, but is a small effect here. Using MKL.jl might also help, instead of the default OpenBLAS library.

Further transposing the data in C, which if I’m thinking correctly ought to bring us closer to what Numpy is actually doing, gives:

```julia
julia> Ct = transpose(permutedims(C)); # same values, row major

julia> Ct == C
true

julia> @btime sum(log, $pa * $Ct);
  210.394 μs (2 allocations: 43.08 KiB)

```

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 26, 2021, 8:59am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/4 "2021-01-26T08:59:29Z")

</div>

Julia is column-major order. Hence you should consider a column vector `p` and a transposed `C` and compute `C * p` in your log-likelihood function:

```julia
using BenchmarkTools
using LinearAlgebra

log_likelihood_no_types(p, C) = sum(log.(p * C))
log_likelihood_no_types_transpose(p, C) = sum(log.(C * p))

Nx = 5500
Np = 200

p = range(0, 1000, length=Np)' |> collect
C = rand(Np, Nx)
pt = collect(p')
Ct = Matrix(C')

# Compile
log_likelihood_no_types(p, C)
log_likelihood_no_types_transpose(pt, Ct)

println("Julia with Nx=$Nx, Np=$Np; BLAS threads: $(BLAS.get_num_threads())")

@btime log_likelihood_no_types($p, $C)
@btime log_likelihood_no_types_transpose($pt, $Ct)

```

I get:

```julia
julia> @btime log_likelihood_no_types($p, $C)
  538.149 μs (4 allocations: 86.16 KiB)
59499.619126650425

julia> @btime log_likelihood_no_types_transpose($pt, $Ct)
  245.567 μs (4 allocations: 86.16 KiB)
59499.619126650425

```

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [January 26, 2021, 9:02am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/5 "2021-01-26T09:02:19Z")

</div>

Comparisons with NumPy seem to trigger immediate replies 🤣.

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [January 26, 2021, 9:03am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/6 "2021-01-26T09:03:22Z")

</div>

Woah, I’m now getting `12.213 μs (2 allocations: 8.12 KiB)` instead of 101! So `adjoint(::Vector{Float64})` is apparently faster than a row-vector… Makes sense, since Julia uses column-major order. But such a drastic difference, holy cow!

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [January 26, 2021, 9:04am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/7 "2021-01-26T09:04:51Z")

</div>

Yep, just changing to `adjoint` and column-major order makes a _huge_ difference, wow!

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 26, 2021, 9:07am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/8 "2021-01-26T09:07:16Z")

</div>

Just to give you a feeling:

```julia
M = rand(1000,1000);

function fcol(M)
    for col in 1:size(M, 2)
        for row in 1:size(M, 1)
            M[row, col] = 42
        end
    end
    nothing
end

function frow(M)
    for row in 1:size(M, 1)
        for col in 1:size(M, 2)
            M[row, col] = 42
        end
    end
    nothing
end

```

gives

```julia
julia> @btime fcol($M)

  673.501 μs (0 allocations: 0 bytes)

julia> @btime frow($M)
  1.422 ms (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [January 26, 2021, 9:07am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/9 "2021-01-26T09:07:20Z")

</div>

Need to show off Julia’s performance against NumPy? Pythonistas hate this one simple trick:

```
p = collect(range(0, 1000, length=Np))'

```

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [January 26, 2021, 9:09am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/10 "2021-01-26T09:09:19Z")

</div>

Duuuuude… Wait, I’m gonna test if NumPy’s performance on column vectors will be as slow as Julia’s with row vectors…

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 26, 2021, 9:10am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/11 "2021-01-26T09:10:40Z")

</div>

Well, if you write the same loops in Python (i.e. do not use numpy) you’re doomed anyways 🙂

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [January 26, 2021, 9:14am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/12 "2021-01-26T09:14:46Z")

</div>

OK, so I changed `p` and `C` to column-major order in NumPy:

```
p = np.array(np.linspace(0, 1000, Np), order='F', copy=True)
C = np.array(np.random.rand(Np, Nx), order='F', copy=True)

```

And even got a speedup:

```
Python with Nx=5500, Np=200
288 µs ± 6.37 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

```

But Julia got 251.563 μs anyway, so cool! Still strange that Julia got such a speedup from simple reordering. Yes, cache locality and ummm contiguity, but NumPy is _always_ fast, while Julia depends on such quirks…

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [January 26, 2021, 9:17am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/13 "2021-01-26T09:17:59Z")

</div>

Note that the timings you compare are not the same. In Julia `@btime` reports the fastest timing, whereas Python seems to report the mean.

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [January 26, 2021, 9:21am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/14 "2021-01-26T09:21:58Z")

</div>

Huh, I was initially using `@benchmark`, but then saw `@btime` somewhere and switched to that without reading the docs.

With `println(mean(@benchmark log_likelihood_no_types($p, $C)))` I’m now getting:

```
Julia with Nx=5500, Np=200; BLAS threads: 4
TrialEstimate(394.292 μs)

```

…which is still _waaaaay_ better than the _best_ time of 1238 μs I had before.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [January 26, 2021, 9:30am UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/15 "2021-01-26T09:30:12Z")

</div>

> [@ForceBru](#):
>
> `C = np.array(np.random.rand(Np, Nx), order='F', copy=True)`

Transposing the data here is I think like `Ct` → `C` in my answer, and for me it slows down Python (and Julia).

I don’t think `order='F'` changes `p` at all, it’s still a (row) vector, so you still call a routine for matrix-vector products. The original big slowdown came from calling a (more complicated) routine for matrix-matrix multiplication, with one row. I can’t seem to make Python make this mistake, I guess it is checking the sizes at run-time, while Julia is going by the types.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 26, 2021, 1:42pm UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/16 "2021-01-26T13:42:58Z")

</div>

It probably would be smart to automatically handle this. It should be pretty cheap.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 26, 2021, 5:51pm UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/17 "2021-01-26T17:51:43Z")

</div>

Does `log` SIMD nowadays or you need LoopVectorization or Simdpirates etc. for that still?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 26, 2021, 6:00pm UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/18 "2021-01-26T18:00:17Z")

</div>

At the very least, you need `@simd ivdep`, and might need `LoopVectorization`. most of our special functions should be relatively simdable, but the problem is that they need to get inlined to simd, and the functions aren’t marked `@inline` because you only want to inline them when you can vectorize them. Furthermore `log` and 64 bit `exp` are table based, you you need `ivdep` for `@simd` to realize that the table accesses don’t break things. `@avx` is probably easier at this point.

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [January 26, 2021, 6:31pm UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/19 "2021-01-26T18:31:29Z")

</div>

I can’t find the paper at the moment, but there was a paper written about `BenchmarkTools` where the authors showed that minimum time is actually a better estimator of performance than the mean. This is from the [docs](https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md):

> The minimum is a robust estimator for the location parameter of the time distribution,  
> and should not be considered an outlier  
> The median, as a robust measure of central tendency, should be relatively unaffected  
> by outliers  
> The mean, as a non-robust measure of central tendency, will usually be positively  
> skewed by outliers  
> The maximum should be considered a primarily noise-driven outlier, and can change  
> drastically between benchmark trials.

So you might want to just consider the minimum runtime for each benchmark.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [January 26, 2021, 11:39pm UTC](https://discourse.julialang.org/t/sum-log-p-c-is-2-to-4-times-slower-than-in-numpy/53960/20 "2021-01-26T23:39:03Z")

</div>

If you’re doing this repeatedly, it pays to preallocate your output vector, and you can get some easy bonus gains with `vmapreduce` from LoopVectorization:

```julia
julia> @btime sum(log, $Ct*$pt)
  199.300 μs (2 allocations: 43.08 KiB)
59500.644576618346

julia> tmp = similar(pt, size(Ct, 1));

julia> @btime vmapreduce(log, +, mul!($tmp, $Ct, $pt))
  148.600 μs (0 allocations: 0 bytes)
59500.64457661832

```