Vectorising Objective Function

Question

I've got two functions which are causing a bottleneck in my program. It's for a graph-based combinatorial optimisation problem which works out the average energy by calculating the energy for each bit string.

At the moment, I'm mapping the function mwis_objective over the entire list of bitstrings. The best way I found elsewhere on stackexchange, I've already used!

The statevector is delivered as a dictionary but then is converted to two numpy arrays. The delay doesn't come from this though, but the mwis_objective function itself, particuarly the two for loops in it.

I think there's probably better ways to use the dictionary as well, but I'm not as sure on this: would passing in the lists of G.edges() and G.nodes once speed things up ?

Unfortunately, it's part of a wider program so I can't use multithreading or things like this.

Here is the code at the moment:

def compute_mwis_energy_sv(counts, G):
    '''
    Computes objective value from an inputted dictionary of amplitudes and bit strings and Graph
    Optimised using function mapping from https://stackoverflow.com/questions/14372613/efficient-creation-of-numpy-arrays-from-list-comprehension-and-in-general/14372746#14372746
    '''


    bit_strings  = list(counts.keys())
    amplitudes = np.array(list(counts.values()))
    number_of_mwis_values = len(bit_strings)
    objective_values =np.fromiter((mwis_objective(bit_string,G) for bit_string in bit_strings),float,number_of_mwis_values)
    probabilities = np.abs(amplitudes)**2

    objective  = np.sum(probabilities * objective_values)

    return objective

def mwis_objective(x,G):
    '''
    Takes in networkx graph  G and a bit string  x from the Qasm output and calculates the < psi | C | psi >
    Need to take note of the order of the bit string.
    '''
    array_of_x = np.array(list(x), dtype=int) # this takes the bit string 1001 to a numpy array for faster access
    # getting the maximum weight of nodes
    node_weights = G.nodes(data='node_weight')  # this is how the node weight is stored in my graph attributes
    just_weights = np.array([weight[1] for weight in node_weights]) #gets just the node weight and converts it to a np array
    scale = np.amax(just_weights) # gets the maximum weight so the node weights are
    scaled_weights = just_weights/scale  # used as J_i,j must be greater than weight of node; all node weights are scaled to below 0 and J_ij is put as 2

    objective = 0
    for i,j in G.edges():  # independent set
        if array_of_x[i] == 1 and  array_of_x[j]==1:  # interconnecting nodes are in the same set

         
            objective += 2


    for i in G.nodes():
        if array_of_x[i] == 1:
            objective -=scaled_weights[i]
            

    return objective

And this is what a typical counts dictionary looks like, as well as a networkx function that generates the graphs I've used for completion:

counts = {'000': (0.3357980114755203-0.3765264103419055j), '001': (-0.45109872283358193+0.08796189074046101j), '010': (0.10490787465961775+0.04222801037227937j), '011': (-0.1929723237399522-0.01534995215062387j), '100': (-0.4510987228335819+0.08796189074046087j), '101': (-0.08891853199461548-0.4236712656325255j), '110': (-0.19297232373995227-0.015349952150623658j), '111': (-0.14362094081740967-0.1650614674350345j)}

def weighted_path_graph(number_of_nodes,graph_weights):
    """
    Creates a weighted path graph of default length three with different weights
    Graph_weights is a list with the desired node weights

    """


    path = nx.path_graph(number_of_nodes)
    graph_weights=[1,3,1]
    for i in range(0,number_of_nodes):
        path.nodes[i]["node_weight"] = graph_weights[i]

    return path

Edit:

I've found a way to vectorise the bottom part of the objective function, but am still stuck on the edge list part.

    weights_to_subtract = array_of_x * scaled_weights
    total_weight = np.sum(weights_to_subtract)
    
    objective = objective - total_weight

Answer 1

# you can convert the edges to an array 
edges = np.array(G.edges())

# and then use fancy indexing to find the matches. 
# summing up the resulting array gives you the nubmer of matches
# and multiplying by two counts each match double.
2* np.sum(((array_of_x[edges[:,0]]==1) & (array_of_x[edges[:,1]]==1)))

edit: more detailed explanation:

# this uses every source node in every edge as an index for your array_of_x 
# and returns an array of the values of array_of_x at all these indices
array_of_x[edges[:,0]]

# this uses every target node in every edge as an index for your array_of_x 
# and returns an array of the values of array_of_x at all these indices
array_of_x[edges[:,1]]


# this checks where the return value is exactly one and returns an array of True/False
array_of_x[edges[:,0]]==1

# this does all of the above and additionally checks whether the values are the same
(array_of_x[edges[:,0]]==1) & (array_of_x[edges[:,1]]==1)


# the above line returns an array of True/False.
# when you call np.sum on this, True is interpreted as 1 and False is interpreted as 0. So summing up gives you the number of matches. 
# since in your code you want to count each match double, you can multiply everything by 2:

2 * np.sum(((array_of_x[edges[:,0]]==1) & (array_of_x[edges[:,1]]==1)))

# how this would give negative values is beyond me.