Confusion about prod()

What am I misunderstanding about prod() that poly2 gives a reducing over empty collection error?

# WORKS
function poly(i, x, xx)
    p = 1
    for j in 1:i-1
        p *= (x-xx[j])
    end
    return p
end

# REDUCING OVER EMPTY COLLECTION ERROR
function poly2(i, x, xx)
    return prod(x-xx[j] for j in 1:i-1)
end

EDIT: I’m not entirely certain how to appropriately type annotate the function signatures, but i is a positive integer, x is a lambda, and xx is a list of Float64

Can you give an example of inputs that produce the error? Also telling us what you expect it to return would be helpful.

You say "x is a lambda". What do you mean by that?

I edited just now to explain the function signatures. Re: lambdas, I guess they’re just called anonymous functions, I was thinking of Python.

What I’m doing here is passing in an anonymous function and looping over a list of x values in order to build a polynomial. Specifically these are Newton’s Divided Difference polynomials, and poly is just the helper function.

Here’s the relevant code:

# Newton's divided difference polynomial
function ndd(; data)
    n = length(data)
    # initialize nxn matrix
    m = zeros(Float64, n, n)
    
    # unzips list of (x,y) coords into separate lists
    # set first column of DD matrix to f(x) for x in xx
    xx, m[:,1] = map(y -> map(x -> x[y], data), 1:2)

    # generate lower triangular NDD matrix
    for i in 2:n
        for j in 2:i
            m[i,j] = (m[i,j-1]-m[i-1,j-1])/(xx[i]-xx[i-j+1])
        end
    end
    
    # return an anonymous function in x: 
    # sum the diagonal entries as coefficients
    # on the generated Newton polynomial
    return x-> sum(m[i,i]*poly(i, x, xx) for i in 1:n)
end

It needs a lot of work, as there’s a lot of silly things happening, but this was a past homework assignment that I’m just trying to clean up.

What is the product of an empty collection? prod() will try to return something sensible (for example, for a vector of numbers it will give the multiplicative identity, i.e. 1). But for a generator like x-xx[j] for j in 1:i-1 there may not be a correct answer in the case that 1:i-1 is an empty range.

Your poly() function works because it handles the empty case by returning 1. You can achieve the same thing with the slightly more verbose but equivalent foldl function:

julia> f(x) = foldl(*, x, init=1)
f (generic function with 1 method)

julia> f(1:5)
120

julia> f(1:1)
1

julia> f(2:1) # 2:1 is an empty range
1
1 Like

But note also that there’s nothing wrong with just writing out the for loop. Loops are fast in Julia, so if they fit your problem well then feel free to use them.

1 Like

Ah, okay yeah the first iteration is from 1:0. I was using it to short circuit out of the for loop so that I wouldn’t have to write explicit code for the first term of the divided difference polynomial.

I’d appreciate some pointers on how to make the posted code more “Julian”, if anybody feels up to it.

This is basically the same as

The issue is that when you have an empty product (or an empty sum), the correct answer is “1” (or “0” for sums), where “1” is the multiplicative identity for the operands (“0” is the additive identity for sums). However, it needs to return the multiplicative identity of the correct type. For example, if you are taking a product of integers, the multiplicative identity is 1, but if if you are taking the product of 2x2 matrices then the identity is the 2x2 identity matrix. Julia has a function one(T) to compute the multiplicative identity for a given type T.

If you give the prod function an empty collection as an argument, it tries to use the eltype (element type) of the collection as the argument to one, in order to find the right type of multiplicative identity. Unfortunately, generator expressions (foo for x in X) don’t have an element type that can be inferred by the compiler (https://github.com/JuliaLang/julia/issues/18695).

One solution is to use something like foldl or mapreduce that allows you to explicitly pass the “initial” value of the reduction, which is also the return value in the empty case.

3 Likes

Thank you, this is very helpful. I hope its okay that I marked the previous answer by @rdeits as the solution.