# Numpy 10x faster than Julia ?! What am I doing wrong ?! \[solved - julia faster now\]
**URL:** https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922
**Category:** Performance
**Tags:** question
**Created:** [October 15, 2019, 4:07am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922 "2019-10-15T04:07:37Z")
**Posts on this page:** 18
**Page:** 2
### Author:  [@improbable22](https://discourse.julialang.org/u/improbable22)
#### Post date: [October 15, 2019, 9:19am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/22 "2019-10-15T09:19:05Z")
> [@Sukera](#):
>
> Doing the “matmul” and then the mul with `[1.,0.]` amounts to just the upper left corner needing to be calculated.
My python is rusty, but I thought `r = m1*m0*np.array([1., 0.])` was scalar multiplication, identical to `r = m1 .* m0 .* [1,0]` in Julia. It calculates a 2x2 array `r`, and surely it does all the calculations, even the ones which give zero and are discarded. (And of course scalar `.*` gives the same results as matrix multiplication `*` on diagonal matrices.)
It’s a little hard to tell from this example which of these calculations you actually need done in the real case. If it’s just `m1[1,1] * m0[1,1]` then just do that, don’t make matrices.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [October 15, 2019, 9:21am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/23 "2019-10-15T09:21:23Z")
> [@joak](#):
>
> Exact, this is not the real computation, just a silly example going through the same steps.
Some morals to keep in mind:
1. If you find yourself creating and operating on zillions of tiny arrays in a loop in Julia, consider using StaticArrays
2. If you find yourself creating zillions of big arrays in a loop in Julia, consider working in-place.
3. If you find yourself writing “vector” style code that applies simple functions (e.g. your `f(p)`) one by one to a sequence of big arrays in Julia, consider rewriting in a different style. A “scalar” style that does lots of operations on a _single_ input is fast in Julia (unlike Python), and may be faster than “vector” style code because of memory locality and similar considerations.
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [October 15, 2019, 9:23am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/24 "2019-10-15T09:23:43Z")
yes, it’s scalar, but since the elements at those positions (0,0) and (1,1) are arrays and not a proper materialized diagonal matrix (not even in type), the upper left arrays of m0 and m1 end up being multiplied - I tried this on my laptop, but I’m currently on my phone so I can only show you later. Basically m1[0][0]\*m0[0][0] is the same as the multiplication of t1 and t0 after the lambdas from above are applied.
That doesn’t mean that numpy isn’t calculating it though, just that because of the later r[0] everything else is discarded anyway and is thus not needed.
---
### Author:  [@joak](https://discourse.julialang.org/u/joak)
#### Post date: [October 15, 2019, 9:33am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/25 "2019-10-15T09:33:48Z")
Woaw !
This is my first contact with julia community.
I love it.
So, problem solved: I have now something faster than numpy ( thanks to @davidbp )
And probably even faster with StaticArrays (haven’t try it yet) thanks to @stevengj
- other improvements (thanks to all)
And even more: a bug detected in my python numpy code
I was trying julia, now I am moving to julia
Thanks !
---
### Author:  [@xiaodai](https://discourse.julialang.org/u/xiaodai)
#### Post date: [October 15, 2019, 9:35am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/26 "2019-10-15T09:35:06Z")
> [@joak](#):
>
> I have now something faster than numpy
Maybe ask in the numpy community and see if there can be improvement for a fair comparision?
---
### Author:  [@improbable22](https://discourse.julialang.org/u/improbable22)
#### Post date: [October 15, 2019, 9:37am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/27 "2019-10-15T09:37:51Z")
Yes that’s what I said, mostly. Except that if the question is about speed, then calculating one multiplication, vs. scalar multiplication of two arrays, vs matrix multiplication of two arrays, these are very different things! Even if what you keep from them is the same number.
This is why having a sufficiently nontrivial example that the answer depends on such details would be useful. I still think the Python example skips matrix multiplication entirely, so it’s not a fair comparison:
```julia
>>> np.array([[1,2],[3,4]]) * np.array([[5,6],[7,8]])
array([[5, 12],
[21, 32]])
```
```julia
julia> [1 2; 3 4] .* [5 6; 7 8] # N^2 operations
2×2 Array{Int64,2}:
5 12
21 32
julia> [1 2; 3 4] * [5 6; 7 8] # N^3 operations
2×2 Array{Int64,2}:
19 22
43 50
```
> [@Sukera](#):
>
> the upper left arrays of m0 and m1 end up being multiplied
You mean the elements, surely. m0 is a 2x2 array of complex numbers.
---
### Author:  [@joak](https://discourse.julialang.org/u/joak)
#### Post date: [October 15, 2019, 9:39am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/28 "2019-10-15T09:39:06Z")
> [@xiaodai](#):
>
> > [@joak](#):
> >
> > I have now something faster than numpy
>
> Maybe ask in the numpy community and see if there can be improvement for a fair comparision?
There are many other reasons to move to julia. Same order of magnitude than numpy is enough to choose julia for the project.
---
### Author:  [@xiaodai](https://discourse.julialang.org/u/xiaodai)
#### Post date: [October 15, 2019, 9:40am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/29 "2019-10-15T09:40:11Z")
> [@joak](#):
>
> Same order of magnitude than numpy is enough to choose julia for the project.
How do you know that numpy can not be even faster?
---
### Author:  [@joak](https://discourse.julialang.org/u/joak)
#### Post date: [October 15, 2019, 9:46am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/30 "2019-10-15T09:46:58Z")
[/quote]
How do you know that numpy can not be even faster?
[/quote]
Maybe numpy can be slightly faster for some operations but not for all, because:
1. julia has jit
2. julia’s ability to handle multicore
3. (and gpu)
---
### Author:  [@xiaodai](https://discourse.julialang.org/u/xiaodai)
#### Post date: [October 15, 2019, 9:49am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/31 "2019-10-15T09:49:25Z")
> [@joak](#):
>
> - julia has jit
> - julia’s ability to handle multicore
> - (and gpu)
Numba
---
### Author:  [@joak](https://discourse.julialang.org/u/joak)
#### Post date: [October 15, 2019, 9:57am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/32 "2019-10-15T09:57:01Z")
> [@xiaodai](#):
>
> > [@joak](#):
> >
> > - julia has jit
> > - julia’s ability to handle multicore
> > - (and gpu)
>
> Numba
Numba is very limited.
Pypy is great if you want python and do not need numpy.
Python3 is good for notebooks (a huge collections of libraries available).
For heavy computations Julia seems to be the best.
Confirmed here.
---
### Author:  [@davidbp](https://discourse.julialang.org/u/davidbp)
#### Post date: [October 15, 2019, 10:04am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/33 "2019-10-15T10:04:41Z")
I would not state confirmed here. I actually don’t want to start another “Why not numba/cython/shedskin/pypy/tensorflow” argument in this thread.
You can make python with Cython fast but… you have to rewrite your code, type everything (event loop indicies) or the performance will degrade a lot. If you are in a situation where you **only care about speeding up functions that get as input numpy arrays** then cython is fine. Once you start working with dataframes, text, python lists that contain different types… then it’s another whole new story. You will have to learn C++ and how to use C++ data structures in cython. And obviously… compile your code and link external libraries etc…
If Python users end up in this discussion I will leave a link to some
[cython\_experiments](https://github.com/davidbp/python_tutorials/blob/master/python_advanced/cython/02_cython_numpy.ipynb).
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [October 15, 2019, 10:05am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/34 "2019-10-15T10:05:35Z")
> [@improbable22](#):
>
> You mean the elements, surely. m0 is a 2x2 array of complex numbers.
Yes the elements, but they’re arrays!
```python
>>> n= 10**6
>>> p = 2 * np.pi * np.random.rand(2,n)
>>> t0,t1 = p
>>> type(t0)
>>> m0 = np.array([[cos(t0) - 1j*sin(t0), 0], [0, cos(t0) + 1j*sin(t0)]])
>>> type(m0)
>>> m0
array([[array([ 0.89155991+0.45290277j, -0.89750168-0.44101105j,
-0.90581293+0.42367786j, ..., 0.50766908+0.86155215j,
0.93237733+0.36148652j, 0.97745514-0.2111432j ]),
0],
[0,
array([ 0.89155991-0.45290277j, -0.89750168+0.44101105j,
-0.90581293-0.42367786j, ..., 0.50766908-0.86155215j,
0.93237733-0.36148652j, 0.97745514+0.2111432j ])]], dtype=object)
```
Notice the dtype being generic `object`, indicating that m0 holds a mixture of types.
And thus:
```python
>>> m1 = np.array([[cos(t1) - 1j*sin(t1), 0], [0, cos(t1) + 1j*sin(t1)]])
>>> def z(x):
... return cos(x) - 1*j*sin(x)
...
>>> (m1[0][0]*m0[0][0]==z(t1)*z(t0)).all()
True
>>> ((m1*m0)[0][0]==z(t1)*z(t0)).all()
True
```
To be clear here, I agree with you! It’s scalar multiplication all the way down, I’m just saying that because of that and the fact that the later `r[0]` basically ignores everthing but the upper left cell, there’s room for not even doing any other muls than those needed for the upper left cell.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 15, 2019, 10:46am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/35 "2019-10-15T10:46:04Z")
> [@joak](#):
>
> I am writing a quantum computer simulator, for something like this [https://qml.entropicalabs.com/](https://qml.entropicalabs.com/) (this one is in js)
You might want to take a look at [GitHub - QuantumBFS/Yao.jl: Extensible, Efficient Quantum Algorithm Design for Humans.](https://github.com/QuantumBFS/Yao.jl)
And I would definitely use static arrays here as Steven suggested.
---
### Author:  [@improbable22](https://discourse.julialang.org/u/improbable22)
#### Post date: [October 15, 2019, 11:21am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/36 "2019-10-15T11:21:33Z")
> [@Sukera](#):
>
> Yes the elements, but they’re arrays!
OK, my bad, I mis-read the Python.
But I still think there is no matrix multiplication at all on the Python side here. Unless I’m still misreading things? This is confusing since everybody is comparing to Julia with matrix multiplications.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 15, 2019, 11:21am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/37 "2019-10-15T11:21:41Z")
There is Pythran too. It seems to have threaded loop fusion which we dont have (yet)…
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [October 15, 2019, 11:48am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/38 "2019-10-15T11:48:22Z")
Yes you’re right, there’s no matmuls or even vecmuls here, just scalar multiplication. I’ve misread that in my first post up top as well, hence the misplaced `mul!` there.
---
### Author:  [@xiaodai](https://discourse.julialang.org/u/xiaodai)
#### Post date: [October 15, 2019, 11:52am UTC](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922/39 "2019-10-15T11:52:25Z")
> [@rveltz](#):
>
> threaded loop fusion
What is that?
[Previous page](https://discourse.julialang.org/t/numpy-10x-faster-than-julia-what-am-i-doing-wrong-solved-julia-faster-now/29922.md?page=1)