how to vectorize cvxpy code to speed up compilation?

Question

I'm using cvxpy on Linux with mosek-10.0.37 solver.

I'm running the code below and I'm having a issue when I increase the value of K (after 40). The cvxpy gives the warnings below and takes a long time in the reduction and consumes a lot of memory before crashing. Please, could you give some advice to solve this problem?

cvxpy UserWarning: Constraint #2 (and #3) contains too many subexpressions. Consider vectorizing your CVXPY code to speed up compilation.

Regards, Jose

========

import numpy as np
import cvxpy as cp

K = 5

e_gp = {(0, 0): 0.01, (0, 1): 0.01, (0, 2): 0.01, (0, 3): 0.01, (0, 4): 0.01, (1, 0): 0.01, (1, 1): 0.01, (1, 2): 0.01, (1, 3): 0.01, (1, 4): 0.01, (2, 0): 0.01, (2, 1): 0.01, (2, 2): 0.01, (2, 3): 0.01, (2, 4): 0.01, (3, 0): 0.01, (3, 1): 0.01, (3, 2): 0.01, (3, 3): 0.01, (3, 4): 0.01, (4, 0): 0.01, (4, 1): 0.01, (4, 2): 0.01, (4, 3): 0.01, (4, 4): 0.01}
f_gp = {(0, 0): 0.0022083736662195513, (0, 1): 1.6361300004264737e-07, (0, 2): 6.673065375868695e-07, (0, 3): 1.8216969961107378e-05, (0, 4): 7.747569454751729e-06, (1, 0): 0.0007991844189121604, (1, 1): 0.0958477227950784, (1, 2): 0.0006366120217946077, (1, 3): 0.005370938447086471, (1, 4): 0.0007190951574825933, (2, 0): 0.0001461922258284424, (2, 1): 8.781670998829123e-06, (2, 2): 0.013922894884250001, (2, 3): 0.0007882566113873586, (2, 4): 5.609072687532707e-05, (3, 0): 1.6647180495599578e-05, (3, 1): 9.701212848634301e-06, (3, 2): 8.415947875397424e-05, (3, 3): 0.002340059870467834, (3, 4): 6.260839608528006e-06, (4, 0): 0.002396438491055845, (4, 1): 2.777751959882801e-05, (4, 2): 0.00011551904109760868, (4, 3): 0.0005849270607107239, (4, 4): 0.013415428308115004}
r_gp = {0: 3.9230406757174874e-05, 1: 0.04336223814165744, 2: 0.0009650237242044663, 3: 5.605495010604486e-05, 4: 0.0008930969228883886}

sinr_t = cp.Variable(pos=True)

p = cp.Variable(shape=(K,), pos=True)
p_min = 1e-10 * np.ones(K)
p_max = 1.0 * np.ones(K)

t_e, t_f, t_r = [], [], []
for k in range(K):
    t_e.append(cp.sum(cp.hstack(e_gp[(k,kp)]*p[kp] for kp in range(K) if k != kp)))
    t_f.append(cp.sum(cp.hstack(f_gp[(k,kp)]*p[kp] for kp in range(K))))
    t_r.append(r_gp[k])
   
objective = cp.Maximize(sinr_t)
S = (cp.hstack(t_e)+cp.hstack(t_f)+cp.hstack(t_r))/p
constraints = [p >= p_min, p <= p_max, S <= 1/sinr_t]
problem = cp.Problem(objective, constraints)
sol = problem.solve(solver=cp.MOSEK, gp=True, verbose=False)