Fix noisy_simulator panic when initializing with non-square Kraus matrices (microsoft#1826)

When initializing an `Operation`, the first thing we do is transposing
the matrices, `nalgebra` panics internally when transposing non-square
matrices. This PR moves the transposition operation after the dimensions
check to gracefully handle that case and adds a unit test to check for
this in the future.

Notes:
 1. This bug was introduced in the performance improvement PR.
2. It doesn't affect well-formed programs; it just affects programs
using ill-formed Kraus operators.

Author: orpuente-MS

Diff for: noisy_simulator/src/operation.rs

@@ -53,20 +53,6 @@ impl Operation {
     /// Matrices must be of dimension 2^k x 2^k, where k is an integer.
     /// Returns `None` if the kraus matrices are ill formed.
     pub fn new(mut kraus_operators: Vec<SquareMatrix>) -> Result<Self, Error> {
-        // Performance note: Because `nalgebra` stores its matrices in column major
-        // form, we use `gemv_tr` in the `apply_kernel` function when multiplying to avoid
-        // incurring cache misses. That is why we transpose all Kraus operators when they
-        // enter the simulator:
-        //   `gemv(1, matrix, vec, 0)` is equivalent to `gemv_tr(1, matrix_tr, vec, 0)`,
-        // but the later has much better performance.
-        //
-        // SAFETY of transposing: all these matrices are only consumed by the
-        // `kernel.rs/apply_kernel` function which effectively transposes them
-        // back when multypling, so it is safe to do this transformation.
-        for kraus_operator in &mut kraus_operators {
-            kraus_operator.transpose_mut();
-        }
-
         let (dim, _) = kraus_operators
             .first()
             .ok_or(Error::FailedToConstructOperation(
@@ -91,6 +77,20 @@ impl Operation {
             }
         }
 
+        // Performance note: Because `nalgebra` stores its matrices in column major
+        // form, we use `gemv_tr` in the `apply_kernel` function when multiplying to avoid
+        // incurring cache misses. That is why we transpose all Kraus operators when they
+        // enter the simulator:
+        //   `gemv(1, matrix, vec, 0)` is equivalent to `gemv_tr(1, matrix_tr, vec, 0)`,
+        // but the later has much better performance.
+        //
+        // SAFETY of transposing: all these matrices are only consumed by the
+        // `kernel.rs/apply_kernel` function which effectively transposes them
+        // back when multypling, so it is safe to do this transformation.
+        for kraus_operator in &mut kraus_operators {
+            kraus_operator.transpose_mut();
+        }
+
         // Performance note: The effect_matrix = Σᵢ (kᵢ† ⋅ k), but due to performance
         // reasons described above we want to store its transpose. But (A ⋅ B)^T = B^T ⋅ A^T.
         // Therefore, we have to swap the order of the factors.

Diff for: noisy_simulator/src/operation/tests.rs

@@ -46,6 +46,20 @@ fn check_operation_number_of_qubits_is_computed_correctly() {
     assert_eq!(1, op.number_of_qubits());
 }
 
+#[test]
+fn check_non_square_kraus_operator_does_not_panic() {
+    let op = operation!(
+        [
+            0., 0.;
+        ]
+    );
+
+    assert!(matches!(
+        op,
+        Err(crate::Error::FailedToConstructOperation(_))
+    ));
+}
+
 /// Check that the inner matrices of the instrument are constructed correctly.
 #[test]
 fn check_effect_matrix_is_computed_correctly() {