# Kron vs scalar product speed difference. python code faster?

**URL:** <https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [January 13, 2017, 12:24pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460 "2017-01-13T12:24:27Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [January 13, 2017, 2:59pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/21 "2017-01-13T14:59:57Z")

</div>

Final result on `ajf/rowvector/af9a28f`

```julia
 This RBM used scalar product
        BenchmarkTools.Trial: 
  memory estimate: 1.42 mb
  allocs estimate: 814
  --------------
  minimum time: 72.793 ms (0.00% GC)
  median time: 73.184 ms (0.00% GC)
  mean time: 74.461 ms (0.05% GC)
  maximum time: 83.571 ms (0.00% GC)
  --------------
  samples: 68
  evals/sample: 1
  time tolerance: 5.00%
  memory tolerance: 1.00%

```

Edit: Same results for with or without interpolation into `@benchmark`.

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [January 13, 2017, 3:02pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/22 "2017-01-13T15:02:25Z")

</div>

Thank you for taking the time and making an effort to explain things to me. I learned a few new things today, which makes this a good day in my book.

> [@stevengj](#):
>
> The reasons is that the parameters are often (as here) global variables, and you don’t want to benchmark the slow type inference on the global variables. Interpolating basically means that the parameters are evaluated before benchmarking, so that once the benchmark starts Julia already knows their types.

I think I understand this aspect well enough,… or maybe not. The point I don’t get is how globals benefit one and hurt the other.

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [January 13, 2017, 3:14pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/23 "2017-01-13T15:14:36Z")

</div>

> [@mkborregaard](#):
>
> It would be super cool if you could post the final optimized code (to be compared to the code in the OP)

[https://gist.github.com/Evizero/49dbc204b772ed63f6fbbf573e8a21da](https://gist.github.com/Evizero/49dbc204b772ed63f6fbbf573e8a21da)

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [January 13, 2017, 3:17pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/24 "2017-01-13T15:17:25Z")

</div>

Thanks!  
 ← 20 char word limit →

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [January 13, 2017, 3:58pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/25 "2017-01-13T15:58:18Z")

</div>

As a partial aisde, I’ve been working on some code lately in which this naive implementation of `sigmoid` (which would be called `invlogit` in the statistics world) is exactly the main source of problems: this exact definition uses a large proportion of total time in the inner loop of my code and it’s also the primary source of numeric issues because of the limited range of `x` for which it generates a result that’s not exactly 0 or 1.

Also worth noting that StatsFuns.jl implements this function as `logistic` using exactly this naive formulation with some added generic typing: [https://github.com/JuliaStats/StatsFuns.jl/blob/ed4867460e4df2cb3ffd52dcde3875b4b635572f/src/basicfuns.jl#L12](https://github.com/JuliaStats/StatsFuns.jl/blob/ed4867460e4df2cb3ffd52dcde3875b4b635572f/src/basicfuns.jl#L12)

In statistics, you could probably do even better than optimizing `sigmoid` by noting that this function is essentially always mixed with other functions: you typically end up needing a mixture of `log(sigmoid(x))` and `log(1 - sigmoid(x))` for fitting logistic regression models, so optimizing those compositions could provide even greater improvements.

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [January 13, 2017, 6:49pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/26 "2017-01-13T18:49:38Z")

</div>

Thank you for taking the time to look at the code. I didn’t know about the @view option or about the fact that .+= avoids realocating memory. I have tested @Evizero code but I still get a lot of GC. This is what I get

## memory estimate: 1.33 gb allocs estimate: 22845

minimum time: 601.147 ms (20.87% GC)  
median time: 633.693 ms (20.14% GC)  
mean time: 624.172 ms (20.33% GC)  
maximum time: 648.480 ms (19.94% GC)

Should I test the code in julia 0.6? I am testing it in julia 0.5.

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [January 13, 2017, 6:56pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/27 "2017-01-13T18:56:10Z")

</div>

> [@davidbp](#):
>
> memory estimate: 1.33 gb

I get the same result on the last version in 0.5. (the last version is really 0.6 focused)

for 0.6 master: `memory estimate: 3.86 mb`

---

<div class="post-metadata">

**Author:** ![dfdx](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dfdx/32/120_2.png) [@dfdx](https://discourse.julialang.org/u/dfdx)\
**Post date:** [January 13, 2017, 8:38pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/28 "2017-01-13T20:38:14Z")

</div>

> [@stevengj](#):
>
> should really be:
> 
> Delta\_W .+= lr .\* (ehp .\* x’ .- ehn .\* xneg’)

I’m surprised nobody has mentioned BLAS yet - I’ve implemented exactly the [same model](https://github.com/dfdx/Boltzmann.jl/blob/master/src/rbm.jl) a while ago, and BLAS gave pretty good performance boost and zero memory allocation (at least at these earlier days of Julia). For example, in real-life settings the expression above would work with mini-batches of, say, 1000 vectors (e.g. `size(ehp) == (255, 1000); size(x) == (784, 1000)`) and can be rewritten as:

```julia
@benchmark Delta_W = lr .* (ehp * x' .- ehn * xneg')     
 BenchmarkTools.Trial:       
   memory estimate: 5.38 mb 
   allocs estimate: 22   
   --------------                                                
   minimum time: 12.316 ms (0.00% GC)    
   median time: 29.460 ms (0.00% GC)      
   mean time: 28.002 ms (0.76% GC)       
   maximum time: 50.396 ms (0.00% GC)   
   --------------                  
   samples: 179    
   evals/sample: 1     
   time tolerance: 5.00%    
   memory tolerance: 1.00%

```

or using BLAS:

```julia
@benchmark begin                                        
     gemm!('N', 'T', 1.0, ehp, x, 0.0, buf)              
     gemm!('N', 'T', 1.0, ehn, xneg, 0.0, Delta_W)  
     axpy!(-1.0, Delta_W, buf)        
     scal!(176400, lr, Delta_W, 1)   
 end                                          
 BenchmarkTools.Trial:              
   memory estimate: 0.00 bytes 
   allocs estimate: 0 
   --------------                                                
   minimum time: 10.181 ms (0.00% GC)
   median time: 10.623 ms (0.00% GC)  
   mean time: 11.759 ms (0.00% GC)    
   maximum time: 45.104 ms (0.00% GC)     
   --------------                         
   samples: 426           
   evals/sample: 1             
   time tolerance: 5.00%     
   memory tolerance: 1.00%

```

Although I haven’t checked it on 0.6, so latest Julia may still beat these results.

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [January 13, 2017, 8:56pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/29 "2017-01-13T20:56:47Z")

</div>

EDIT: Made a mistake which resulted in unfair comparison. The gist I initially posted here in this exact post was wrong as the BLAS implementation is intended to do the whole loop at once, while the other do just one iteration. This mistake does not affect anything I posted above

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [January 13, 2017, 9:24pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/30 "2017-01-13T21:24:43Z")

</div>

Ah. I think I have a conceptual error here and am comparing the wrong things

EDIT: this is only in reference to my last post concerning BLAS

---

<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:** [January 13, 2017, 10:46pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/31 "2017-01-13T22:46:30Z")

</div>

@Evizero, it might be nice to have some version of this as a benchmark in the BaseBenchmarks package, to help use track Julia performance, if you want to do a pull request.

---

<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:** [January 13, 2017, 10:53pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/32 "2017-01-13T22:53:53Z")

</div>

> [@dfdx](#):
>
> For example, in real-life settings the expression above would work with mini-batches of, say, 1000 vectors (e.g. size(ehp) == (255, 1000); size(x) == (784, 1000)) and can be rewritten as:

As soon as you are working with multiple vectors at once, so that you are doing BLAS3 operations (matrix multiplications etc.), then I agree that you definitely want to exploit a fast BLAS (like the OpenBLAS that Julia links). However, you don’t generally need to call low-level BLAS functions like `gemm!` directly. `A_mul_Bt!` and similar high-level functions are just as efficient.

For BLAS1 operations like `axpy!` and `scal!`, you are probably better off with the fusing broadcast operations. The increase in locality and reduction of other overheads that you get from fusing the loops will beat the minor optimizations that OpenBLAS can do for BLAS1 operations.

(And if most of your time is spent in BLAS3 operations, then you should expect essentially the same performance in Julia and NumPy, assuming that they are using the same BLAS library.)

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [January 14, 2017, 9:39pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/33 "2017-01-14T21:39:20Z")

</div>

Thank you a lot for your time. It turns out that most of the speed came from a transpose that your code avoided, not from the fact that you used A\_mul\_B!..

In your version you wrote  
Delta\_W .+= lr .\* (ehp .\* x’ .- ehn .\* xneg’)

In my version I had  
Delta\_W .+= lr \* ( x \* ehp’ - xneg \* ehn’)’

Here there are some tests  
[https://github.com/davidbp/learn\_julia/blob/master/speed\_tests/comparing\_functions.ipynb](https://github.com/davidbp/learn_julia/blob/master/speed_tests/comparing_functions.ipynb)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 14, 2017, 9:45pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/34 "2017-01-14T21:45:17Z")

</div>

> [@davidbp](#):
>
> Thank you a lot for your time. It turns out that most of the speed came from a transpose that your code avoided, not from the fact that you used A\_mul\_B!..

I think you’d be happy to hear that Julia in v0.6 will “take transposes seriously” via the type system, making it essentially a free operation. It was a huge thread, culminating in:

> <https://github.com/JuliaLang/julia/issues/4774#issuecomment-269142361>
>
> from @alanedelman:
> 
> We really should think carefully about how the transpose of …a vector should dispatch the various \`A\_\*op\*\_B\*\` methods. It must be possible to avoid new types and ugly mathematics. For example, vector'vector yielding a vector (#2472, #2936), vector' yielding a matrix, and vector'' yielding a matrix (#2686) are all bad mathematics.
> 
> What works for me mathematically (which avoids introducing a new type) is that for a 1-dimensional \`Vector\` \`v\`:
> \- \`v'\` is a no-op (i.e. just returns \`v\`),
> \- \`v'v\` or \`v'\*v\` is a scalar,
> \- \`v\*v'\` is a matrix, and
> \- \`v'A\` or \`v'\*A\` (where \`A\` is an \`AbstractMatrix\`) is a vector
> 
> A general \_N\_-dimensional transpose reverses the order of indices. A vector, having one index, should be invariant under transposition. 
> 
> In practice \`v'\` is rarely used in isolation, and is usually encountered in matrix-vector products and matrix-matrix products. A common example would be to construct bilinear forms \`v'A\*w\` and quadratic forms \`v'A\*v\` which are used in conjugate gradients, Rayleigh quotients, etc.
> 
> The only reason to introduce a new \`Transpose{Vector}\` type would be to represent the difference between contravariant and covariant vectors, and I don't find this compelling enough.

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [January 14, 2017, 9:49pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/35 "2017-01-14T21:49:46Z")

</div>

I already tried the code in 0.6 and the penalty of the transpose seemed huge.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 14, 2017, 9:50pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/36 "2017-01-14T21:50:36Z")

</div>

Less than 24 hours ago?

---

<div class="post-metadata">

**Author:** ![davidbp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidbp/32/463_2.png) [@davidbp](https://discourse.julialang.org/u/davidbp)\
**Post date:** [January 14, 2017, 9:52pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/37 "2017-01-14T21:52:01Z")

</div>

Downloaded 0.6 dev 3 hours ago… (macOS)

Julia Version 0.6.0-dev.2069  
Commit ff9a949 (2017-01-13 02:17 UTC)  
Platform Info:  
OS: macOS (x86\_64-apple-darwin13.4.0)  
CPU: Intel(R) Core™ i7-4650U CPU @ 1.70GHz  
WORD\_SIZE: 64  
BLAS: libopenblas (USE64BITINT DYNAMIC\_ARCH NO\_AFFINITY Haswell)  
LAPACK: libopenblas64\_  
LIBM: libopenlibm  
LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 14, 2017, 9:54pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/38 "2017-01-14T21:54:04Z")

</div>

But how old is the nightly you used?

[https://status.julialang.org/](https://status.julialang.org/)

Most nightlies haven’t had an update since then (only COPR has), so if you didn’t build it from source or use COPR, you won’t have that update. It will say how many days from the last master at the top of the REPL when you open it up.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 14, 2017, 9:55pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/39 "2017-01-14T21:55:32Z")

</div>

> [@davidbp](#):
>
> Julia Version 0.6.0-dev.2069

That’s a day old master and won’t have it. I’d retry when it’s in the nightly.

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [January 14, 2017, 10:07pm UTC](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460/40 "2017-01-14T22:07:46Z")

</div>

@davidbp is talking about the outer matrix transpose. I don’t think the rowvector PR affects that

EDIT: also, it looks like `lr * (...)'` gets translated into `A_mul_Bc(lr, ...)`, so I don’t think (?) that the outer transpose influences much in your particular code.

> [@davidbp](#):
>
> In your version you wrote  
> `Delta_W .+= lr .* (ehp .* x' .- ehn .* xneg')`
> 
> In my version I had  
> `Delta_W .+= lr * ( x * ehp' - xneg * ehn')'`

Aside from that, also observe how your version uses `*` and the other `.*`. The `*` (without the dot) causes broadcast fusion to stop, which will cause the creation of temporary arrays. So really there are a few non-obvious sources of performance penalties.

The change @ChrisRackauckas is talking about affects the inner vector transposes and so gives an additional performance boost.

[Previous page](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460.md?page=1)

[Next page](https://discourse.julialang.org/t/kron-vs-scalar-product-speed-difference-python-code-faster/1460.md?page=3)