The problem is that the big-M method works well for ILPs while the variable, here the $p$s, are non-integer. When g_i(h) is very close to 1, I do not get the correct value for h.
I tried it out with different precision levels of CVX but it did not help either.

I was wondering if there is a trick for this in CVX.

I’m not sure whether you are trying to model a strict linear inequality constraint, or a conditional (logical) constraint. Either way, you will have to address the matter of solver feasibility tolerance.

CVX treats and sends to the solver all < constraints as <=. In order to get strict inequality, you will need to introduce a small positive number, small_number. Then use the constraint sum(p(i,:),2) + small_number <= 1 .

In order to allow small_number to be as small as possible, while still dealing with solver tolerance, you will need to adjust the solver feasibility tolerance (applicable to linear (inequality) constraints) from its default value. Do this using cvx_solver_settings( '{name}', {value} ) . See http://cvxr.com/cvx/doc/solver.html#advanced-solver-settings . You will need to examine the solver documentation to determine the parameter name to adjust.

I believe that cvx_precision controls solver termination criteria, but not solver feasibility tolerance, which I believe is is what you need to adjust.