Performance issue in python- network creation based on the Euclidean distance with nested for loops is too slow

Question

I want to create a network in which the links are formed based on a similarity metric defined as the Euclidean distance between the nodes. The distance is calculated using socio-demographic features of customers such as gender and age. The problem is the code takes 200 seconds to just create the network and as I am tuning my model and the code executes at least 100 times, the long execution time of this piece is making the whole code run slowly.

So, the nodes are in fact customers. I defined a class for them. They have two attributes gender (numerical; specified by number 0 or 1) and age (varies from 24 to 44) which are stored in a csv file. I have generated a sample csv file here :

#number of customers
ncons = 5000
gender = [random.randint(0, 1) for i in range(ncons)]
age = [random.randint(22, 45) for i in range(ncons)]
customer_df = pd.DataFrame(
    {'customer_gender': gender,
     'customer_age': age
    })
customer_df.to_csv('customer_df.csv', mode = 'w', index=False)

The Euclidean distance delta_ik is of the form following. In the formula, n is the number of attributes (here n=2, age and gender). For customers i and k, S_f,i - S_f,k is the difference between attribute f = 1,2 which is divided by the maximum range of attribute f for all the customers (max d_f). So the distance is the distance in the values of socio-demographic attributes, not geographical positions.

Then I define the similarity metric H_ik which creates a number between 0 and 1 from delta_ik as follow:. Finally, For customers i and k, I generate a random number rho between 0 and 1. If rho is smaller than H_ik, the nodes are connected.

So, the code that keeps delta_ik in a matrix and then uses that to generate the network looks as below:

import random
import pandas as pd
import time
import csv
import networkx as nx
import numpy as np
import math
#Read the csv file containing the part worth utilities of 184 consumers
def readCSVPWU():
    global headers
    global Attr
    Attr = []
    with open('customer_df.csv') as csvfile:
        csvreader = csv.reader(csvfile,delimiter=',')
        headers = next(csvreader)  # skip the first row of the CSV file.
        #CSV header cells are string and should be turned to a float number.
        for i in range(len(headers)):   
            if headers[i].isnumeric():
                headers[i] = float(headers[i])
        for row in csvreader:
            AttrS = row
            Attr.append(AttrS)
    #convert strings to float numbers
    Attr = [[float(j) for j in i] for i in Attr]
    #Return the CSV as a matrix with 17 columns and 184 rows 
    return Attr

#customer class
class Customer:
    def __init__(self, PWU = None, Ut = None):
        self.Ut = Ut
        self.PWU = Attr[random.randint(0,len(Attr)-1)]  # Pick random row from survey utility data  


#Generate a network by connecting nodes based on their similarity metric
def Network_generation(cust_agent):
    start_time = time.time() # track execution time

    #we form links/connections between consumeragentsbasedontheirdegreeofsocio-demographic similarity.
    global ncons
    Gcons = nx.Graph()
    #add nodes
    [Gcons.add_node(i, data = cust_agent[i]) for i in range(ncons)]
    #**********Compute the node to node distance
    #Initialize Deltaik with zero's
    Deltaik = [[0 for xi in range(ncons)] for yi in range(ncons)] 
    #For each attribute, find the maximum range of that attribute; for instance max age diff = max age - min age = 53-32=21
    maxdiff = []
    allval = []
    #the last two columns of Attr keep income and age data
    #Make a 2D numpy array to slice the last 2 columns (#THE ACTUAL CSV FILE HAS MORE THAN 2 COLUMNS)
    np_Attr = np.array(Attr)
    #Take the last two columns, income and age of the participants, respectively
    socio = np_Attr[:, [len(Attr[0])-2, len(Attr[0])-1]]
    #convert numpy array to a list of list
    socio = socio.tolist()
    #Max diff for each attribute

    for f in range(len(socio[0])):
        for node1 in Gcons.nodes():
        #keep all values of an attribute to find the max range
            allval.append((Gcons.nodes[node1]['data'].PWU[-2:][f]))
        maxdiff.append((max(allval)-min(allval)))
        allval = []
# THE SECOND MOST TIME CONSUMING PART ********************

    for node1 in Gcons.nodes():
        for node2 in Gcons.nodes():
            tempdelta = 0
            #for each feature (attribute)
            for f in range(len(socio[0])):
                Deltaik[node1][node2] = (Gcons.nodes[node1]['data'].PWU[-2:][f]-Gcons.nodes[node2]['data'].PWU[-2:][f])
                #max difference
                insidepar = (Deltaik[node1][node2] / maxdiff[f])**2
                tempdelta += insidepar
            Deltaik[node1][node2] = math.sqrt(tempdelta)
     # THE END OF THE SECOND MOST TIME CONSUMING PART ********************
       
    #Find maximum of a matrix
    maxdel = max(map(max, Deltaik))
    #Find the homopholic weight
    import copy
    Hik = copy.deepcopy(Deltaik)
    for i in range(len(Deltaik)):
        for j in range(len(Deltaik[0])):
            
            Hik[i][j] =1 - (Deltaik[i][j]/maxdel)
    #Define a dataframe to save Hik
    dfHik = pd.DataFrame(columns = list(range(ncons) ),index = list(range(ncons) ))
    temp_h = []
    #For every consumer pair $i$ and $k$, a random number $\rho$ from a uniform distribution $U(0,1)$ is drawn and compared with $H_{i,k}$ . The two consumers are connected in the social network if $\rho$ is smaller than $H_{i,k}$~\cite{wolf2015changing}.
# THE MOST TIME CONSUMING PART ********************
    for node1 in Gcons.nodes():
        for node2 in Gcons.nodes():
            #Add Hik to the dataframe
            temp_h.append(Hik[node1][node2])
            rho = np.random.uniform(0,1,1)
            if node1 != node2:
                if rho < Hik[node1][node2]:
                    Gcons.add_edge(node1, node2)
        #Row idd for consumer idd keeps homophily with every other consumer
        dfHik.loc[node1] = temp_h
        temp_h = []
    # nx.draw(Gcons, with_labels=True)            
    print("Simulation time: %.3f seconds" % (time.time() - start_time))
# THE END OF THE MOST TIME CONSUMING PART ********************

    return Gcons     
#%%
#number of customers
ncons = 5000
gender = [random.randint(0, 1) for i in range(ncons)]
age = [random.randint(22, 39) for i in range(ncons)]
customer_df = pd.DataFrame(
    {'customer_gender': gender,
     'customer_age': age
    })
customer_df.to_csv('customer_df.csv', mode = 'w', index=False)
readCSVPWU()
customer_agent = dict(enumerate([Customer(PWU = [], Ut = []) for ij in range(ncons)])) # Ut=[]
G = Network_generation(customer_agent)

I realized that there are two nested loops that are more time consuming than others, but I am not sure how to write them more efficiently. I would tremendously appreciate if you could please give me some advice on the ways to decrease the elapsed time.

Thank you so much

Answer 1

Boy, those are huge numbers. 5000 Nodes results in an adjacency matrix with 25 million entries. But since Deltaik(a,b) == Deltaik(b,a), you don't actually need to consider the whole cartesian product. In fact, your current code doubles the chance for an edge to be created because each pair of nodes (a,b) gets two chances: one for (a,b) and one for (b,a). This would be the first improvement:

Instead of the nested loop:

for node1 in Gcons.nodes:
    for node2 in Gcons.nodes:
        ...

Use itertools.combinations:

for node1, node2 in itertools.combinations(Gcons.nodes, 2):
    ...

Memory consupmtion

There will still be 12497500 combinations, so let's talk about memory consumption. You are creating huge temporary data structures (Deltaik, Hik) that are only iterated once and then discarded. I'd suggest the following datastructure:

delta_iks = []

That's really all you need. The reason is that itertools.combinations is insanely fast, and it will generate the combinations in the exact same order every time (as long as you don't modify the graph, of course). So whenever you need to iterate all edges, i.e. (node1, node2), just generate them again. You also don't have to create a separate list of Hik because you can calculate it on the fly when you actually add the edges to the graph.

Runtime

12497500 combinations also means that the calculation of the distances will be executed that many times, so you should make sure that it is efficient.

First of all, there are a lot of unnecessary dictionary accesses only to retrieve the customer instance of each node. You can access them more directly. Also, the [-2:] is unnecessary because PWU already has only 2 elements, but it generates a new sublist each time it is called.

for (node1, customer1), (node2, customer2) in itertools.combinations(Gcons.nodes.data('data'), 2):
...
    customer1.PWU[f]...

for f in range(len(socio[0])):

This is a red flag in python. What you really want to do in that loop is simultaneously iterate the features of customer1 and customer2 and the maxdiff of that feature. Since all of them are stored in the same order, this can be achieved using zip:

for c1_attr, c2_attr, maxdiff_attr in zip(customer1.PWU, customer2.PWU, maxdiff):
    tempdelta += ((c1_attr - c2_attr) / maxdiff_attr)**2

The complete first loop now looks like this (and completes in about 7 seconds on my computer):

delta_iks = []
for (_, customer1), (_, customer2) in itertools.combinations(Gcons.nodes.data('data'), 2):
    tempdelta = 0
    #for each feature (attribute)
    for c1_attr, c2_attr, maxdiff_attr in zip(customer1.PWU, customer2.PWU, maxdiff):
        tempdelta += ((c1_attr - c2_attr) / maxdiff_attr)**2
    delta_iks.append(math.sqrt(tempdelta))
max_delta_ik = max(delta_iks)

As you can see, we don't even use the nodes' names anymore. If we need to iterate all edges (node1, node2) with their corresponding delta_ik, we can simply do:

for (node1, node2), delta_ik in zip(itertools.combinations(Gcons.nodes, 2), delta_iks):
    ...

Now prepare the Hik. Note that this is a generator that will only generate each value when you iterate it. It's not a list because that would again use a lot of memory.

h_iks = (1-(delta_ik / max_delta_ik) for delta_ik in delta_iks)

One thing about numpy.random.uniform that I also didn't know before: generating one value at a time 12 million times is significantly slower (20 seconds on my machine) than telling numpy how many values you are going to need beforehand. So let's prepare the rhos now:

rhos = np.random.uniform(0,1,len(delta_iks))

networx.Graph has a very convenient method for adding multiple edges: add_edges_from, which takes any iterable, e.g. a generator. So let's create yet another generator that yields all edges if a random rho is less than its h_ik:

edges = (edge for edge, h_ik, rho in zip(itertools.combinations(Gcons.nodes, 2), h_iks, rhos) if rho < h_ik)

Notice how nothing has been calculated yet, and the only significant amount of memory is used by the delta_iks list. Now add the edges to the graph (this takes about 8-9 seconds on my machine):

Gcons.add_edges_from(edges)

The complete Network_generation function:

#Generate a network by connecting nodes based on their similarity metric
def Network_generation(cust_agent):
    start_time = time.time() # track execution time

    #we form links/connections between consumeragentsbasedontheirdegreeofsocio-demographic similarity.
    global ncons
    Gcons = nx.Graph()
    #add nodes
    [Gcons.add_node(i, data = cust_agent[i]) for i in range(ncons)]
    #**********Compute the node to node distance
    #For each attribute, find the maximum range of that attribute; for instance max age diff = max age - min age = 53-32=21
    maxdiff = []
    allval = []
    #the last two columns of Attr keep income and age data
    #Make a 2D numpy array to slice the last 2 columns (#THE ACTUAL CSV FILE HAS MORE THAN 2 COLUMNS)
    np_Attr = np.array(Attr)
    #Take the last two columns, income and age of the participants, respectively
    socio = np_Attr[:, [len(Attr[0])-2, len(Attr[0])-1]]
    #convert numpy array to a list of list
    socio = socio.tolist()
    #Max diff for each attribute
    for f in range(len(socio[0])):
        for node1 in Gcons.nodes():
        #keep all values of an attribute to find the max range
            allval.append((Gcons.nodes[node1]['data'].PWU[-2:][f]))
        maxdiff.append((max(allval)-min(allval)))
        allval = []

    # THE SECOND MOST TIME CONSUMING PART ********************
    delta_iks = []
    for (_, customer1), (_, customer2) in itertools.combinations(Gcons.nodes.data('data'), 2):
        tempdelta = 0
        #for each feature (attribute)
        for c1_attr, c2_attr, maxdiff_attr in zip(customer1.PWU, customer2.PWU, maxdiff):
            tempdelta += ((c1_attr - c2_attr) / maxdiff_attr)**2
        delta_iks.append(math.sqrt(tempdelta))
    max_delta_ik = max(delta_iks)
     # THE END OF THE SECOND MOST TIME CONSUMING PART ********************

    h_iks = (1-(delta_ik / max_delta_ik) for delta_ik in delta_iks)
    rhos = np.random.uniform(0,1,len(delta_iks))
    edges = (edge for edge, h_ik, rho in zip(itertools.combinations(Gcons.nodes, 2), h_iks, rhos) if rho < h_ik)

    # THE MOST TIME CONSUMING PART ********************
    Gcons.add_edges_from(edges)
    # THE END OF THE MOST TIME CONSUMING PART ********************
    print("Simulation time: %.3f seconds" % (time.time() - start_time))
    return Gcons

This takes around 16 seconds on my computer. For reference, your original code took roughly 105 seconds on my computer.

I haven't touched the part before "THE SECOND MOST TIME CONSUMING PART" because it doesn't significantly contribute to memory consumption and run time. It could still be improved, though.

---

Further improvements

Memory

The delta_iks list takes up around 500MB of heap space. When an edge has been added to the graph, we don't need the corresponding delta_ik anymore, so we could theoratically remove it from the list. Unfortunately, python lists don't provide an efficient way to do so. But collections.deques do.

Make delta_iks a deque:

# delta_iks = []
delta_iks = collections.deque()

Instead of defining h_iks with a generator expression, we'll make it a generator function that yields the same values, but additionally removes them from delta_iks. Note that now it has to be invoked with () in the edges generator:

# h_iks = (1-(delta_ik / max_delta_ik) for delta_ik in delta_iks)
def h_iks():
    while delta_iks:
        yield(1-(delta_iks.popleft() / max_delta_ik))

rhos = np.random.uniform(0,1,len(delta_iks))
# edges = (edge for edge, h_ik, rho in zip(itertools.combinations(Gcons.nodes, 2), h_iks, rhos) if rho < h_ik)
edges = (edge for edge, h_ik, rho in zip(itertools.combinations(Gcons.nodes, 2), h_iks(), rhos) if rho < h_ik)

The reason why you only save 400MB is that the empty deque is still about 100MB on the heap. Keep in mind, though, that the graph is significantly larger (somewhere in the region of 1.8GB), so it won't make that much of a difference.

Runtime

If I understood your code correctly, the Hik is the probability that an edge is created. If you only need it for that purpose, you could shave off some numerical calculations. Instead of

if rho < (1-(delta_iks / max_delta_ik))

You could test

if rho > (delta_iks / max_delta_ik)

And furthermore, instead of normalizing delta_ik and picking a 0 <= rho < 1, you cold use the "raw" delta_ik and pick a 0 <= rho < max_delta_ik:

rhos = np.random.uniform(0,max_delta_ik,len(delta_iks))
edges = (edge for edge, delta_ik, rho in zip(itertools.combinations(Gcons.nodes, 2), delta_iks, rhos) if rho > delta_ik)

That, however, won't make a significant difference (3 seconds on my machine).