Hacker News Reader

Comment by threatofrain | Hacker News Reader

Optimization will never end up that way unless the system understands both your intention and the domain, but why would you want to express your intentions on matrices through 3 for-loops, as opposed to a product operation already available to you?

And if a matrix product operation were available to a language as a standard construct, wouldn't it be every bit as opaque, encouraging students to think that hand-written for-loops aren't as effective as using black-box things?

View on HN

I am not advocating to use hand-written for-loops for matrix multiplication. Just that the language and the community should not impose an undue penalty to for loops. Not all scientific algoritms are matrix multiplication. Sometimes, the computation you want to perform is easily expressed using loops, and your language should not get in the way. I have seen many students struggling to "vectorize" their algorithm and writing an ugly mess instead of a clean loop. They were under the impression that writing a loop over the whole data was extremely dirty and to be avoided at all costs. When I wrote the thing with a simple loop, it was so ridiculously slow, that it even re-inforced their apprehension to loops. This is a sad tragedy, and there is no real reason for it to be that way.

I've had this experience when writing a simple algorithm for run length encoding, a problem which I struggled to vectorize. Compare the numpy algo with pure python:

  def rle_encode_pure_python(sequence):
    """Run length encode array."""
    previous = sequence[0]
    count = 1
    out = []
    for element in sequence[1:]:
        if element == previous:
            count += 1
        else:
            out.append((count, previous))
            previous = element
            count = 1
    out.append((count, previous))
    return out


  def rle_encode_numpy(sequence):
    diffs = np.concatenate(( np.array((True, )), np.diff(sequence)!=0))
    indices = np.concatenate((np.where(diffs)[0],np.array((sequence.size, ))))
    counts = np.diff(indices).astype('uint16')
    values = sequence[diffs].astype('uint8')
    return np.rec.fromarrays((counts, values),names=('count','value'))

The numpy version is faster, but a lot less legible (and makes some assumptions about data types and array length). Luckily the pure python version can be jitted with numba and is faster than numpy.

This is a great example, I love it, thank you very much!

Note: the numpy version is missing the import.

Note2: I tend to prefer using numba than numpy. Yet, sometimes numpy is inevitable, especially for linear algebra routines. In that case, I am in a world of pain because numba and numpy do not interact well at all!

< numba and numpy do not interact well at all

Can you please elaborate on this? Not challenging you but trying to understand the nuances.

Numba tutorial says otherwise: https://numba.pydata.org/numba-doc/dev/user/5minguide.html "Numba likes NumPy functions"

Nice and concise example. Would love a mini-writeup / benchmarks for this (including the numba version).

Here is the code in a notebook: https://github.com/Python-Meetup-Rotterdam/meetup1/blob/mast...

Gave a little presentation on my findings on using different approaches

Returning a list? Nay, Yield!

Completely agreed, would never write it this way anymore but the example is a bit older

Nested for loops, right? Did you know that going in different order can be much faster as it affects memory access patterns?

That's cool and fun to know, and some people may enjoy writing optimisations like that (and we def need those people) but so many people just want to crunch numbers and make graphs. They should just call the vectorized library functions and get on with their research.

Yeah, Python of today is pretty much what Fortran was for the previous scientific generation.

Reply on news.ycombinator.com