Improve LinearQubitOperator by porting functionality from Cirq

[kevinsung]

Currently, we have a function generate_linear_qubit_operator in utils/_linear_qubit_operator.py which takes a QubitOperator and outputs a LinearQubitOperator or ParallelLinearQubitOperator (subclasses of scipy.sparse.linalg.LinearOperator) which can be used, for instance, in scipy.sparse eigenvalue routines. This class can be made a lot more useful by improving the implementation of LinearQubitOperator._matvec.

The code below uses Cirq and constitutes a functional implementation of LinearQubitOperator._matvec:

def apply_qubit_operator_to_state(qubit_op, state):
    num_qubits = state.shape[0].bit_length() - 1
    result = np.zeros(state.shape, dtype=state.dtype)
    for term, coeff in qubit_op.terms.items():
        mutated_state = np.reshape(np.copy(state), (2,) * num_qubits)
        for index, action in term:
            if action == 'X':
                pauli = cirq.X
            elif action == 'Y':
                pauli = cirq.Y
            else:
                pauli = cirq.Z
            buffer = np.empty(mutated_state.shape, dtype=state.dtype)
            args = cirq.ApplyUnitaryArgs(target_tensor=mutated_state,
                                         available_buffer=buffer,
                                         axes=(index,))
            mutated_state = cirq.apply_unitary(pauli, args)
        result += coeff*np.reshape(mutated_state, state.shape)
    return result

I found that this implementation gives an 80-100x speedup at large (at least 16 qubits) system sizes with smaller but still significant gains at smaller ones.

The task for this issue is basically to port this code over into an implementation of LinearQubitOperator._matvec. The class ParallelLinearQubitOperator can then be deleted since it would be obsolete.