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

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.

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

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.

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

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.

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