
Iterative methods done right (life's too short to write for-loops) - stellalo
http://lostella.github.io/blog/2018/07/25/iterative-methods-done-right
======
wenc
Like the other commenters have remarked, I have a feeling this post is meant
for folks in numerical computation.

In numerical mathematics, iteration is typically expressed as recursive
equations, and I suspect the iterable construct in Julia helps lower the
impedance mismatch of translating numerical ideas to code. There's a lot of
tedious "book-keeping" (keeping track of indices, etc.) involved in numerical
codes, and appropriate abstractions can help reduce errors, as well as improve
composability and extensibility. At the very least, with iterables the code
will match the math better than if for-loops were used.

For most other (general) programmers though, a conjugate-gradient example
might be a little too complicated of an example for demonstrating iterables.
Also, while iterables as an abstraction can be useful at times, the simplicity
of for-loops can't be beat for most general use cases. There is also a
performance overhead with iterables.

Also, just to point out there is a bit of overlapping but differing
terminology here: "iterative methods" [1] refer to a specific procedure in
computational math; whereas "iterables" in computer programming refer to
iterators [2] which perform (lazy or eager) traversal. The latter can be used
to implement the former, but they are mostly different ideas.

[1]
[https://en.wikipedia.org/wiki/Iterative_method](https://en.wikipedia.org/wiki/Iterative_method)

[2]
[https://en.wikipedia.org/wiki/Iterator](https://en.wikipedia.org/wiki/Iterator)

~~~
bunderbunder
> There is also a performance overhead with iterables

Though, for the sake of pedantry, it's not something that can't be handled by
an optimizing compiler that understands them. Haskell, for example, is good at
this.

~~~
celrod
For what it's worth, just sticking to the Fib example, since that's a fairly
standard benchmark:

    
    
      struct FibonacciIterable{I}
        s0::I
        s1::I
      end
      import Base: iterate
      iterate(iter::FibonacciIterable) = iter.s0, (iter.s0, iter.s1)
      iterate(iter::FibonacciIterable, state) = state[2], (state[2], sum(state))
      function iteratorfib(N)
        i = 0
        x = 0
        for F in FibonacciIterable(0,1)
            i += 1
            x = F
            i > N && break
        end
        x
      end
    

vs

    
    
      function loopfib(N)
        x0 = 0; x1 = 1
        N == 0 && return 0
        for i in 2:N
            x0, x1 = x1, x0+x1
        end
        x1
      end
    

Then I get:

    
    
      julia> N = 30;
    
      julia> using BenchmarkTools
    
      julia> @benchmark iteratorfib($N)
      BenchmarkTools.Trial: 
        memory estimate:  0 bytes
        allocs estimate:  0
        --------------
        minimum time:     31.351 ns (0.00% GC)
        median time:      32.723 ns (0.00% GC)
        mean time:        32.891 ns (0.00% GC)
        maximum time:     54.525 ns (0.00% GC)
        --------------
        samples:          10000
        evals/sample:     994
    
      julia> @benchmark loopfib($N)
      BenchmarkTools.Trial: 
        memory estimate:  0 bytes
        allocs estimate:  0
        --------------
        minimum time:     37.831 ns (0.00% GC)
        median time:      37.871 ns (0.00% GC)
        mean time:        38.616 ns (0.00% GC)
        maximum time:     105.626 ns (0.00% GC)
        --------------
        samples:          10000
        evals/sample:     992
    

If instead of writing "N=30" and interpolating it ($N) I just used N, I'd be
benchmarking the dynamic dispatch. If I declared N constant, or just wrote 30
directly, the compiler turns into a cheater:

    
    
      julia> @benchmark iteratorfib(30)
      BenchmarkTools.Trial: 
        memory estimate:  0 bytes
        allocs estimate:  0
        --------------
        minimum time:     2.990 ns (0.00% GC)
        median time:      2.997 ns (0.00% GC)
        mean time:        3.030 ns (0.00% GC)
        maximum time:     56.348 ns (0.00% GC)
        --------------
        samples:          10000
        evals/sample:     1000
    
      julia> @benchmark loopfib(30)
      BenchmarkTools.Trial: 
        memory estimate:  0 bytes
        allocs estimate:  0
        --------------
        minimum time:     2.990 ns (0.00% GC)
        median time:      2.997 ns (0.00% GC)
        mean time:        3.006 ns (0.00% GC)
        maximum time:     8.908 ns (0.00% GC)
        --------------
        samples:          10000
        evals/sample:     1000
    
    

Anyway, my points here are: 1) In this benchmark that should be focused on
iteration, the custom iteration protocol actually seems faster than the for
loop for some reason. 2) While for iterative algorithms like conjugate
gradient, it does seem like a nice fit, it's probably not the clearest
approach everywhere.

~~~
bunderbunder
And Julia, too, apparently.

I've really got to take a closer look at it. How does it compare to Python,
ecosystem-wise? Pandas might be the thing I'd be most likely to miss in the
stuff I'm doing.

~~~
ChrisRackauckas
There's a lot of generic numerical analysis libraries which is where it shines
over Python. It matches pretty well in areas like data science and machine
learning. For example, for tables there JuliaDB which is more akin to dask
since it works out-of-core (but it's also able to use built-in online
algorithms), DataFrames.jl which matches R, and Pandas.jl which is a straight
wrapper for Python pandas.

------
joe_the_user
Iteration is great but these extra-constructs seem like a unnecessary
proliferation of syntactic sugar.

Which is to say, writing a for-loop is two easy lines, why worry about that?

Most of the convenience I see in the code of the OP is a language that can
return pairs or N-tuples of values. That's useful for a multitude of things.

If you do operations with vectors that all the same iteration and subscription
operations, then one-off errors aren't very likely. One-off/fence-post errors
are much more likely when you iterate through a sequence assigned values in
more complex or staggered fashion and here an iterator-type isn't going to
save you since you have to figure out the logic involved (which you might
encapsulate with some object or type or which might have idiosyncratic logic
all its own).

~~~
dm3
I think this post is mostly showing off the new iteration protocol in Julia
0.7/1.0, not advocating for converting all the for loops into decorators.

However, the post doesn't make it clear, so it's just a guess.

------
RobLach
I love for loops because in one line I can grasp what the loop intends to
compute and usually leads me towards what I should pay attention to when
debugging.

This all seems to require me to remember how it’s implemented and to think
about what’s happening, which is more work than I’d like to do. Perhaps the
assumptions you’d bring along when using a language like Julia alleviates it a
bit, but still if someone else looks at my code they’ll have to look through
all this fanciness.

~~~
BeetleB
>I love for loops because in one line I can grasp what the loop intends to
compute and usually leads me towards what I should pay attention to when
debugging.

C++ programmer here.

In my team I've actively pushed _not_ to have for loops when something from
the algorithms library will do (transform, copy_if, find_if, etc). The reason:
Because I usually _cannot_ quickly grasp what the for loop intends to compute
- I'm quite surprised you're saying the opposite.

I met with some resistance, as people aren't used to thinking in terms of
algorithms. So I was slowly browsing our code base, finding for loops, and
converting the trivial ones to calls to the algorithms library. Then I came
across one for loop that was, frankly, a pointless loop. It didn't seem to do
anything useful. I was sure it was a bug, but I couldn't tell from the context
what the code was supposed to do. I emailed the committer, and he fixed it.
Then I took the fixed version and converted it to an algorithms call.

Then I made a presentation to the team showing some of my examples, including
the one with the bug. That bug convinced many to start using other constructs.
Every time you write a for loop for which there is an algorithm in the
library, you are reinventing the wheel, and have a good chance of introducing
a bug.

(Sorry, probably a bit of a tangent to the whole thread).

~~~
srean
Very true, wish your comment had more upvotes to bubble up to the top of the
thread.

Its a lot easier to miss edge cases, have off by one errors when one is
manipulating iter count directly. Working with the abstraction natural to the
problem reduces these. So manipulating /creating / merging rows, columns,
blocks when working with arrays help avoid many of these problems. Usually the
algorithm pseudocode is written in terms of such abstractions.

------
slx26
>> I’ll illustrate this by examples in Julia, using the conjugate gradient
method for positive (semi)definite linear systems as guinea pig

Which links to a paper that says "An Introduction to the Conjugate Gradient
Method _Without the Agonizing Pain_ ". I mean, I understand that's the kind of
things you use Julia for, but as a dumb programmer, I see the article title,
want to read, but I'm a bit thrown off by that kind of things.

I scrolled a bit through the article and at the beginning I had the same
impression other commenters had here. I just see a few loops like in any
modern lanaguage that accepts iterables [for element in iterable], and I'm
quite confused until I reach the end.

Not really dramatic, but I just wanted to mention it, since the poster seems
to be the author, I feel the title, introduction to the topic, point of the
article, examples, etc, could be made clearer/improved for the general
audience. Or maybe I'm just the wrong audience. Anyway.

~~~
stilley2
I think the intended audience is optimization people, but the title didn't
make that clear. Side note: a lot of folks really love that conjugate gradient
reference, but it's not my favorite.

~~~
benrbray
What is your favorite?

------
newen
It probably got missed in the who conjugate gradient mess, but I think the
main point of the post was that instead of using one giant for loop that
contains your logic, logging, stopping conditions, etc., you can seperate out
your for loop into different parts so that you can mix, match, and reuse your
logic code, loggging code, stopping conditions code, etc.

------
sametmax
No knowing julia, I'm not sure what I'm looking at.

Is that the julia equivalent of 'yield' (as seen in c#, python, js, etc.)
being newly implemented and showed off ?

~~~
stabbles
No, it's really just an iterator in Julia 0.7 or 1.0. At the same time it's
just a function call.

    
    
        for x in itr
          println(x)
        end
    

Is exactly equivalent to and syntactic sugar for

    
    
        y = iterate(itr)
        while y !== nothing
          x, state = y
          println(x)
          y = iterate(itr, state)
        end
    

What makes this interesting is two things:

1\. It uses Julia's multiple dispatch: every iterator has to implement or
extend the function `iterate(itr::MyIteratorType, state...) = ...`

2\. Julia can handle small union types efficiently. In the example above the
inferred type of `y` will be `Union{Nothing,Tuple{...,...}}`. The compiler
will generate efficient code -- in previous versions of Julia any `Union` type
would result in a heap allocation and `y` would internally be a pointer. In
Julia 0.7 there's clever code generation and no allocations.

------
still_grokking
One of the first examples, namely:

    
    
      import Base: iterate
      iterate(iter::FibonacciIterable) = iter.s0, (iter.s0, iter.s1)
      iterate(iter::FibonacciIterable, state) = state[2], (state[2], sum(state))
    

looks a bit like Prolog I think. Interesting.

