Combination of matrix elements giving non-zero value (PYTHON)

Question

I have to evaluate the following expression, given two quite large matrices A,B and a very complicated function F: The mathematical expression

I was thinking if there is an efficient way in order to first find those indices i,j that will give a non-zero element after the multiplication of the matrices, so that I avoid the quite slow 'for loops'.

Current working code

# Starting with 4 random matrices 
A = np.random.randint(0,2,size=(50,50))
B = np.random.randint(0,2,size=(50,50))
C = np.random.randint(0,2,size=(50,50))
D = np.random.randint(0,2,size=(50,50))
indices []
for i in range(A.shape[0]):
    for j in range(A.shape[0]):
        if A[i,j] != 0:
            for k in range(B.shape[1]):
                if B[j,k] != 0:
                for l in range(C.shape[1]):
                    if A[i,j]*B[j,k]*C[k,l]*D[l,i]!=0:
                        indices.append((i,j,k,l))
print indices

As you can see, in order to get the indices I need I have to use nested loops (= huge computational time).

care to share your existing code? Please have a look at [mcve]. — jpp, Feb 08 '18 at 16:03
You say "multiplication of the matrices" but your expression doesn't look like matrix multiplication... Are you sure your indices are all correct? — Stefan Pochmann, Feb 08 '18 at 16:14

AGN Gazer · Answer 1 · 2018-02-08T22:10:34.127

My guess would be NO: you cannot avoid the for-loops. In order to find all the indices ij you need to loop through all the elements which defeats the purpose of this check. Therefore, you should go ahead and use simple array elementwise multiplication and dot product in numpy - it should be quite fast with for loops taken care by numpy.

However, if you plan on using a Python loop then the answer is YES, you can avoid them by using numpy, using the following pseudo-code (=hand-waving):

i, j = np.indices((N, M)) # CAREFUL: you may need to swap i<->j or N<->M
fs = F(i, j, z) # array of values of function F
                # for a given z over the index grid
R = np.dot(A*fs, B) # summation over j
# return R # if necessary do a summation over i: np.sum(R, axis=...)

If the issue is that computing fs = F(i, j, z) is a very slow operation, then you will have to identify elements of A that are zero using two loops built-in into numpy (so they are quite fast):

good = np.nonzero(A) # hidden double loop (for 2D data)
fs = np.zeros_like(A)
fs[good] = F(i[good], j[good], z) # compute F only where A != 0

Combination of matrix elements giving non-zero value (PYTHON)

1 Answers1