Iterative methods done right (life's too short to write for-loops) | Hacker News Reader

Iterative methods done right (life's too short to write for-loops)

lostella.github.io

111 pointsstellalo27 comments

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

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

> 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.

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

  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.

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.

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.

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).

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.

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.

I did laugh at my own ignorance a bit when they say "given the simplicity of the method", apparently referring to an incomprehensible (to me) bunch of maths jargon.

I sometimes wish I had better maths, but then it's so rarely I ever need it, I think I'm better off learning another language.

It's not always obvious that you needed it when you don't understand it.

This. The more math you learn, the more relevant it becomes in everything you do. If you can't do calculus, you'll likely never see where it was you could have been using it to your advantage.

this is probably true. Thanks for the motivation to go learn some more :)

I don't think you need to be good at maths to use Julia. I'm glad I learned some basic trigonometry at school, though...

>when they say "given the simplicity of the method"

I think they meant the implementation, in the end, is not complicated: there are a few variables that get updated according to some formulas, and that's it; there's no branching, no need for fancy data structures, etc.

The method itself is not simple - neither deriving it nor fine-tuning is intuitive, as can be said about other variations on the theme of gradient descent.

The theme is simple: to reach a bottom of a hill, go in the steepest downhill direction[1]. Making it work - and making it work fast - requires modifications to this idea, which require quite a bit of work.

[1]https://en.wikipedia.org/wiki/Gradient_descent

In more conventional scripting languages than Julia, I wish more people would write simple for-loops for simple problems, instead of cramming in some functional construct, just because the language allows it. The extra indirection, the added object allocation, the deepening of the call-stack, it all just crufts up the whole routine. For-loops are simple, execute tidily, and keep code localized, that's good stuff.

>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).

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.

Compelling arguments here. What you’re saying is true in my experience as well in the cases you’ve explained.

> Because I usually cannot quickly grasp what the for loop intends to compute - I'm quite surprised you're saying the opposite.

Yeah, a for loop breaks the relevant information up into several lines and depending on language it adds quite a bit of boilerplate.

INitial state, precondition, postcondition and increment right in front of your eyes. For-loops are rather rich in context information

>> 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.

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.

What is your favorite?

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.

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 ?

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.

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.

Reply on news.ycombinator.com