A Functional way to write map and then reduce?

New to Julia
1

I am trying to implement GDQ method in Julia. Writing code in an imperative schema is a straightforward practice but I feel stuck in rewriting it with functional idioms below is the part of code which I am hanged in. Is it possible to rewrite it in functional schema in an efficient way?

Imperative code:

function coeff_1(points)
  m1_ = m1(points)
  a_ij(i, j) = m1_[i]/(points[i] - points[j])/m1_[j]

  count = length(points)
  result = zeros((count, count))

  for i = 1:count
    for j = 1:count
      if (i != j)
        result[i, j] = a_ij(i, j)
        result[i, i] -= result[i, j]
      end
    end
  end
  result
end

Incomplete functional functional implementation:

function coeff_1_functional(points)
  m1_ = m1(points)
  a_ij(i, j) = m1_[i]/(points[i] - points[j])/m1_[j]

  count = length(points)
  
  result = [a_ij(i, j) for i = 1:count, j = 1:count]
  # the below part is where I stuck, It is both imperative and inefficient
  for i = 1:count
     result[i, i] = -sum(result[i, :])
  end
  result
end

With Regards

2

Or maybe this is what you are looking for.

result = [ i!=j?a_ij(i,j):2*a_ij(i,j) for i = 1:count, j = 1:count ] - diagm( [sum([a_ij(i,j) for j =1:count ]) for i = 1:count] )
3

This is just a straightforward expression and may not help with the efficiency issue. I won’t really comment on the efficiency issue.

function coeff_f(points)
    m = m1(points)
    ix = 1:length(points)

    a(i,j) = m[i]/(points[i] - points[j])/m[j]
    d(i) = -sum(a(i,j) for j = ix if i != j)

    [ ((i == j) ? d(i) : a(i,j)) for i = ix, j = ix ]
end

But might it help to define the auxiliary function for the diagonal entries? Put the condition only there and then compute the diagonal separately.

4

By putting the diagonal condition in the auxiliary, I mean this: instead of NaN, make a(i,j) zero on the diagonal, so it can be used more freely.

function coeff_g(points)
    m = m1(points)
    ix = 1:length(points)

    a(i,j) = (i == j) ? 0.0 : m[i]/(points[i] - points[j])/m[j]
    d(i) = -sum(a(i,j) for j = ix)

    [ a(i,j) for i = ix, j = ix ] + diagm([ d(i) for i = ix ])
end
5

Out of curiosity, what is the problem with the first version of the code?

6

If you refer to my two versions, I think the second version is just a nicer expression. Making a(i,j) defined, as 0.0, on the diagonal made it usable without the conditions.

7

You can make transposed versions of m1_ and points, then broadcast:

function set_diagonal!(matrix, vector)
   set_diagonal_at_index!(i) = 
     matrix[i, i] = vector[i]
   set_diagonal_at_index!.(1:length(vector) )
   matrix
end

function coeff_1_functional(points)
  m1_ = m1(points)
  result = m1 ./ ( points .- points') / m1'
  column_sums = sum(a, 2)
  set_diagonal!(result, column_sums)
end

I haven’t tested it out, but I think it’s about right. Best in 0.6 which has lazy transposes

8

I think there is something with diagonal elements.
Maybe the else part of ternary operator should be zero.
Thanks for your time.

9

It’s completely correct and efficient, just for a simple(functional) and efficient implementation.

10

Nice implementation full of tricks.
Thank you.

11

Yes, that’s so nice.
I just tested it, but it seems to be memory hungry.
Is it the trade of being functional, or there is a better way?

12

The following should also do the work. It is both short and hopefully quick. Unfortunately, because m_1 isn’t detailed in the question, it is hard to make sure everything is 100% correct. So, I’ve googled GDQ to look for the value of m_1 and added a function to calculate it which might be cleaner than the current implementation you have.

This also defines the useful setdiag! function, but makes it both mutating and returning the mutated object, which plays nicely in composition with further processing.

The use of inv.(m1) is important because doing 1/x is slow on most computers and is better done only once per element of m1.

The code:

setdiag!(m,v) = begin foreach(i->m[i,i]=v[i],1:length(v)) ; m end

m_1(points) = prod(setdiag!(points.-points',ones(length(points))),2)

function coeff_1_functional(points)
  m1 = m_1(points)
  result = setdiag!(m1.*inv.(points.-points').*inv.(m1)',zeros(length(points)))
  return setdiag!(result,-sum(result,2))
end
13

I had implemented m_1 as this:

function m1(points)
    function m1_map(row_point)
        row_reduce(prod_, point) = ifelse(point != row_point, (row_point - point) * prod_, prod_)
        reduce(row_reduce, 1, points)
    end

    map(m1_map, points)
end

The most interesting thing about your vectorized implementation is that parts are putted together in a nice way, some code reuse but again it is not as efficient as imperative one. Maybe here the best way is to use C like code.

Thank you.

14

Now with the m1 implementation, I could check both functions give the same results (Yep!). Benchmarking shows my m_1 is 10x faster. As for the coeff_1, the loop code in the original question 2x faster than my coeff_1_functional in Julia 0.5 and about the same on Julia 0.6 (which enjoys better optimization of vectorized operations).

PS benchmark used points vector 10 elements long.