ORTools Linear Prog with Sparsity

from ortools.linear_solver import pywraplp
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

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

diff = np.linalg.inv(A[S_bar,:][:,S_bar] - nu * np.eye(len(S_bar)))
prod = np.matmul(diff, A[S_bar,:][:,S])

# The above is $prod = (A_{\bar{S}, \bar{S}} - \nu I_{|S|})^{-1} A_{\bar{S}, S}$

solver = pywraplp.Solver.CreateSolver('GLOP')
u_S_prog = np.array([solver.NumVar(0, solver.infinity(), "u_S_" + str(i)) for i in range(len(S))])
z = np.array([solver.NumVar(0, solver.infinity(), 'z_' + str(i)) 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)
for index, z_i in enumerate(z):
    solver.Add(z_i <= rhs[index]) # Each component of z must be ≤ each component of rhs

objective = solver.Objective()
for z_i in z:
    objective.SetCoefficient(z_i, 1)

objective.SetMaximization()
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print('Objective value =', solver.Objective().Value())
    print()
    print('Problem solved in %f milliseconds' % solver.wall_time())
    print('Problem solved in %d iterations' % solver.iterations())
    print('Problem solved in %d branch-and-bound nodes' % solver.nodes())
else:
    print(f'The problem with status {status} does not have an optimal solution.')