PULP (GLPK) Linear Prog + Sparsity

import pulp as pl
import pulp
import numpy as np
import scipy.sparse

A = scipy.sparse.random(100, 100, density=0.1).todense()
indices = list(range(100))
S = np.random.permutation(indices)[:10]
S_bar = np.array(list(set(indices) - set(S)))
nu = 0.55
u_S_budget = 50 # Add some scalar here

# Zero out rows of seed nodes
A[S[:, None], S] = 0
A[S[:, None], S_bar] = 0

# $prod = (A_{\bar{S}, \bar{S}} - \nu I_{|S|})^{-1} A_{\bar{S}, S}$
diff = np.linalg.inv(A[S_bar,:][:,S_bar] - nu * np.eye(len(S_bar)))
prod = np.matmul(diff, A[S_bar,:][:,S])

solver = pl.getSolver('GLPK_CMD')
prob = pulp.LpProblem('objective', pulp.LpMaximize)

# Both variables range from 0 to infinity.
u_S_prog = np.array([pl.LpVariable("u_S_" + str(i), 0) for i in range(len(S))])
z = np.array([pl.LpVariable('z_' + str(i), 0) for i in range(len(S_bar))])

rhs = -prod.dot(u_S_prog) # rhs of our to-be-formed minimum linear constraint
rhs = np.array(rhs).squeeze(0)

#prob += sum(x for x in u_S_prog)
prob += sum(x for x in z)

# Add the budget constraints on u_S_prog: needed for non-zero feasible and optimal solution to be possible.
prob += (sum(x for x in u_S_prog) <= u_S_budget)

for index, z_i in enumerate(z):
    prob += (z_i <= rhs[index])

status = prob.solve(solver)
print(prob.objective.value())