faster SLSQP in measurement mitigation (code included)

[aeddins-ibm]

What is the expected behavior?

CompleteMeasFitter.filter.apply was a bottleneck for my qiskit project, so I tried optimizing it.

apply calls Scipy's minimize using SLSQP to infer what the counts "should have been" given the known measurement error statistics. The cost function provided to minimize does some matrix multiplication. This is pretty fast, but given the large number of parameters (one for each state label, which I think will be ~2^N or thereabouts), I suspect (?) the gradient must be expensive (~2^N evaluations). Minimize gives an option to provide the gradient (jacobian) directly, so if the preceding logic is valid then it should speed up the fit routine by ~2^N. This is roughly consistent with a small test I ran for 5 qubits where it decreased the runtime for correcting some hundreds of circuits from ~5-10 seconds to <1 second. More testing would be a good idea (and maybe double-checking my pen-and-paper derivative is right...).

Also added the jacobian for the constraint, not sure what effect this has but seemed easy enough.

I also changed the initial guess, formerly it was a list of small random values, now it's set to the measured (uncorrected) counts, since hopefully the errors will usually be small-ish so that should be a pretty good guess for the corrected values.

Code is below. Sorry it's annoying that this is an issue and not a PR, but as a qiskit user as opposed to a qiskit dev, I wanted to contribute this optimization but don't have the bandwidth right now to learn how to submit a PR, how to install the master branch, etc. Open to suggestions for low-barrier ways for users to give feedback.

New code:

qiskit/ignis/mitigation/measurement/filters.py

raw_data_row = raw_data2[data_idx]
nshots = sum(raw_data_row)
cal_mat = self._cal_matrix
nlabels = len(raw_data_row)
def fun(estimated_corrected_data):
        return np.sum( (raw_data_row - cal_mat.dot(estimated_corrected_data) )**2)

def gradient(estimated_corrected_data):
        return( 2*(cal_mat.dot(estimated_corrected_data) - raw_data_row).dot(cal_mat))
            
cons = ({'type': 'eq', 'fun': lambda x: nshots - np.sum(x), 'jac': lambda x: -1*np.ones_like(x)})
bnds = tuple((0, nshots) for x in raw_data_row)
res = minimize(fun, raw_data_row, method='SLSQP', constraints=cons, bounds=bnds, tol=1e-6, jac=gradient)

Old code:

nshots = sum(raw_data2[data_idx])
def fun(x):
        return sum((raw_data2[data_idx] - np.dot(self._cal_matrix, x))**2)
x0 = np.random.rand(len(self._state_labels))
x0 = x0 / sum(x0)
cons = ({'type': 'eq', 'fun': lambda x: nshots - sum(x)})
bnds = tuple((0, nshots) for x in x0)
res = minimize(fun, x0, method='SLSQP', constraints=cons, bounds=bnds, tol=1e-6)
raw_data2[data_idx] = res.x