Product of a sequence in NumPy

I need to implement this following function with NumPy -

where F_l(x) are N number of arrays that I need to calculate, which are dependent on an array G(x), that I am given, and A_j are N coefficients that are also given. I would like to implement it in NumPy as I would have to calculate F_l(x) for every iteration of my program. The dummy way to do this is by for loops and ifs:

import numpy as np
A = np.arange(1.,5.,1)
G = np.array([[1.,2.],[3.,4.]])

def calcF(G,A):
    N = A.size
    print A
    print N
    F = []
    for l in range(N):
        F.append(G/A[l])
        print F[l]
        for j in range(N):
            if j != l:
                F[l]*=((G - A[l])/(G + A[j]))*((A[l] - A[j])/(A[l] + A[j]))
    return F

F= calcF(G,A)
print F

As for loops and if statements are relatively slow, I am looking for a NumPy witty way to do the same thing. Anyone has an idea?

Listed in this post is a vectorized solution making heavy usage of NumPy's powerful broadcasting feature after extending dimensions of input arrays to 3D and 4D cases with np.newaxis/None at various places according to the computation involved. Here's the implementation -

# Get size of A
N = A.size

# Perform "(G - A[l])/(G + A[j]))" in a vectorized manner
p1 = (G - A[:,None,None,None])/(G + A[:,None,None])

# Perform "((A[l] - A[j])/(A[l] + A[j]))" in a vectorized manner
p2 = ((A[:,None] - A)/(A[:,None] + A))

# Elementwise multiplications between the previously calculated parts
p3 = p1*p2[...,None,None]

# Set the escaped portion "j != l" output as "G/A[l]"
p3[np.eye(N,dtype=bool)] = G/A[:,None,None]
Fout = p3.prod(1)

# If you need separate arrays just like in the question, split it
Fout_split = np.array_split(Fout,N)

Sample run -

In [284]: # Original inputs
     ...: A = np.arange(1.,5.,1)
     ...: G = np.array([[1.,2.],[3.,4.]])
     ...: 

In [285]: calcF(G,A)
Out[285]: 
[array([[-0.        , -0.00166667],
        [-0.01142857, -0.03214286]]), array([[-0.00027778,  0.        ],
        [ 0.00019841,  0.00126984]]), array([[  1.26984127e-03,   1.32275132e-04],
        [ -0.00000000e+00,  -7.93650794e-05]]), array([[-0.00803571, -0.00190476],
        [-0.00017857,  0.        ]])]

In [286]: vectorized_calcF(G,A) # Posted solution 
Out[286]: 
[array([[[-0.        , -0.00166667],
         [-0.01142857, -0.03214286]]]), array([[[-0.00027778,  0.        ],
         [ 0.00019841,  0.00126984]]]), array([[[  1.26984127e-03,   1.32275132e-04],
         [ -0.00000000e+00,  -7.93650794e-05]]]), array([[[-0.00803571, -0.00190476],
         [-0.00017857,  0.        ]]])]

Runtime test -

In [289]: # Larger inputs
     ...: A = np.random.randint(1,500,(400))
     ...: G = np.random.randint(1,400,(20,20))
     ...: 

In [290]: %timeit calcF(G,A)
1 loops, best of 3: 4.46 s per loop

In [291]: %timeit vectorized_calcF(G,A)  # Posted solution 
1 loops, best of 3: 1.87 s per loop

Vectorization with NumPy/MATLAB : General approach

Felt like I could throw in my two cents on my general approach and I would think others follow similar strategies when trying to vectorize codes, especially in a high level platform like NumPy or MATLAB. So, here's a quick check-list of things that could be considered for Vectorization -

Idea about extending the dimensions : Dimensions are to be extended for the input arrays such that the new dimensions hold results that would have otherwise gotten generated iteratively within the nested loops.

Where to start vectorizing from? Start from the deepest (that loop where the code is iterating the most) stage of computation and see how inputs could be extended and the relevant computation could be brought in. Take good care of tracing the iterators involved and extend dimensions accordingly. Move outwards onto outer loops, until you are satisfied with the vectorization done.

How to take care of conditional statements? For simple cases, brute force compute everything and see how the IF/ELSE parts could be taken care of later on. This would be highly context specific.

Are there dependencies? If so, see if the dependencies could be traced and implemented accordingly. This could form another topic for discussion, but here are few examples I got myself involved with.