Add ApplyUnitary operation

Standard/src/Synthesis/MatrixUtils.cs

@@ -0,0 +1,78 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+using System.Diagnostics;
+using System.Numerics;
+using Microsoft.Quantum.Simulation.Core;
+
+namespace Microsoft.Quantum.Synthesis
+{
+    internal class MatrixUtils
+    {
+        /// <summary>
+        /// Checks whether given matrix is unitary.
+        /// </summary>
+        public static bool IsMatrixUnitary(Complex[,] matrix, double tol = 1e-10)
+        {
+            int n = matrix.GetLength(0);
+            if (matrix.GetLength(1) != n) return false; // Unitary matrix must be square.
+
+            for (int i = 0; i < n; i++)
+            {
+                for (int j = 0; j < n; j++)
+                {
+                    Complex dotProduct = 0.0;
+                    for (int k = 0; k < n; k++)
+                    {
+                        dotProduct += matrix[i, k] * Complex.Conjugate(matrix[j, k]);
+                    }
+                    Complex expectedDotProduct = (i == j) ? 1.0 : 0.0;
+                    if ((dotProduct - expectedDotProduct).Magnitude > tol)
+                    {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+
+        /// <summary>
+        /// Converts matrix from C# array to Q# array.
+        /// </summary>
+        public static QArray<QArray<Quantum.Math.Complex>> MatrixToQs(Complex[,] b)
+        {
+            long n1 = b.GetLength(0);
+            long n2 = b.GetLength(1);
+            var a = new QArray<Quantum.Math.Complex>[n1];
+            for (long i = 0; i < n1; i++)
+            {
+                var row = new Quantum.Math.Complex[n2];
+                for (int j = 0; j < n2; j++)
+                {
+                    row[j] = new Quantum.Math.Complex((b[i, j].Real, b[i, j].Imaginary));
+                }
+                a[i] = new QArray<Quantum.Math.Complex>(row);
+            }
+            return new QArray<QArray<Quantum.Math.Complex>>(a);
+        }
+
+        /// <summary>
+        /// Converts matrix from Q# array to C# array.
+        /// </summary>
+        public static Complex[,] MatrixFromQs(IQArray<IQArray<Quantum.Math.Complex>> a)
+        {
+            long n1 = a.Length;
+            long n2 = a[0].Length;
+            var b = new Complex[n1, n2];
+            for (long i = 0; i < n1; i++)
+            {
+                Debug.Assert(a[i].Length == n2);
+                for (long j = 0; j < n2; j++)
+                {
+                    b[i, j] = new Complex(a[i][j].Real, a[i][j].Imag);
+                }
+            }
+            return b;
+        }
+    }
+}

Standard/src/Synthesis/TwoLevelUnitary.cs

@@ -0,0 +1,71 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+using System;
+using System.Diagnostics;
+using System.Numerics;
+using Microsoft.Quantum.Simulation.Core;
+
+namespace Microsoft.Quantum.Synthesis
+{
+    /// <summary>
+    /// Represents a square matrix which is an identity matrix with elements on positions
+    /// (i1, i1), (i1, i2), (i2, i1), (i2, i2) replaced with elements from unitary matrix <c>mx</c>. 
+    /// </summary>
+    internal class TwoLevelUnitary
+    {
+        private Complex[,] mx;   // 2x2 non-trivial principal submatrix.
+        private int i1, i2;      // Indices of non-trivial submatrix.
+
+        public TwoLevelUnitary(Complex[,] mx, int i1, int i2)
+        {
+            Debug.Assert(MatrixUtils.IsMatrixUnitary(mx), "Matrix is not unitary");
+            this.mx = mx;
+            this.i1 = i1;
+            this.i2 = i2;
+        }
+
+        public void ApplyPermutation(int[] perm)
+        {
+            i1 = perm[i1];
+            i2 = perm[i2];
+        }
+
+        // Ensures that index1 < index2.
+        public void OrderIndices()
+        {
+            if (i1 > i2)
+            {
+                (i1, i2) = (i2, i1);
+                (mx[0, 0], mx[1, 1]) = (mx[1, 1], mx[0, 0]);
+                (mx[0, 1], mx[1, 0]) = (mx[1, 0], mx[0, 1]);
+            }
+        }
+
+        // Equivalent to inversion (because matrix is unitary).
+        public void ConjugateTranspose()
+        {
+            mx[0, 0] = Complex.Conjugate(mx[0, 0]);
+            mx[1, 1] = Complex.Conjugate(mx[1, 1]);
+            (mx[0, 1], mx[1, 0]) = (Complex.Conjugate(mx[1, 0]), Complex.Conjugate(mx[0, 1]));
+        }
+
+        // Applies A = A * M, where M is this matrix.
+        public void MultiplyRight(Complex[,] A)
+        {
+            int n = A.GetLength(0);
+            for (int i = 0; i < n; i++)
+            {
+                (A[i, i1], A[i, i2]) = (A[i, i1] * mx[0, 0] + A[i, i2] * mx[1, 0],
+                                        A[i, i1] * mx[0, 1] + A[i, i2] * mx[1, 1]);
+            }
+        }
+
+        public bool IsIdentity(double tol = 1e-10) =>
+            (mx[0, 0] - 1).Magnitude < tol && mx[0, 1].Magnitude < tol &&
+            mx[1, 0].Magnitude < tol && (mx[1, 1] - 1).Magnitude < tol;
+
+        // Converts to a tuple to be passed to Q#.
+        public (IQArray<IQArray<Quantum.Math.Complex>>, long, long) ToQsharp() =>
+            (MatrixUtils.MatrixToQs(this.mx), this.i1, this.i2);
+    }
+}

Standard/src/Synthesis/UnitaryDecomposition.cs

@@ -0,0 +1,170 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+using System;
+using System.Collections.Generic;
+using System.Diagnostics;
+using System.Linq;
+using System.Numerics;
+using Microsoft.Quantum.Simulation.Core;
+
+namespace Microsoft.Quantum.Synthesis
+{
+    internal partial class TwoLevelDecomposition
+    {
+        public class Native : TwoLevelDecomposition
+        {
+            private const double tol = 1e-10;
+
+            public Native(IOperationFactory m) : base(m) { }
+
+            // Returns special unitary 2x2 matrix U, such that [a, b] U = [c, 0],
+            // where c is real and positive. 
+            private static Complex[,] MakeEliminatingMatrix(Complex a, Complex b)
+            {
+                Debug.Assert(a.Magnitude >= tol && b.Magnitude >= tol);
+                double theta = System.Math.Atan((b / a).Magnitude);
+                double lmbda = -a.Phase;
+                double mu = System.Math.PI + b.Phase - a.Phase - lmbda;
+                return new Complex[,] {{
+                    System.Math.Cos(theta) * Complex.FromPolarCoordinates(1.0, lmbda),
+                    System.Math.Sin(theta) * Complex.FromPolarCoordinates(1.0, mu)
+                }, {
+                    -System.Math.Sin(theta) * Complex.FromPolarCoordinates(1.0, -mu),
+                    System.Math.Cos(theta) * Complex.FromPolarCoordinates(1.0, -lmbda)
+                }};
+            }
+
+            // Returns list of two-level unitary matrices, which multiply to A.
+            //
+            // Matrices are listed in application order.
+            // Every matrix has indices differring in 1.
+            // A is modified as result of this function.
+            private static List<TwoLevelUnitary> TwoLevelDecompose(Complex[,] A)
+            {
+                int n = A.GetLength(0);
+                var result = new List<TwoLevelUnitary>();
+
+                for (int i = 0; i < n - 2; i++)
+                {
+                    for (int j = n - 1; j > i; j--)
+                    {
+                        // At this step we will multiply A (from the right) by such two-level 
+                        // unitary that A[i, j] becomes zero.
+                        var a = A[i, j - 1];
+                        var b = A[i, j];
+                        Complex[,] mx;
+                        if (A[i, j].Magnitude <= tol)
+                        {
+                            // Element is already zero. We shouldn't do anything, which equivalent 
+                            // to multiplying by identity matrix.
+                            mx = new Complex[,] { { 1.0, 0.0 }, { 0.0, 1.0 } };
+                            // But if it's last in a row, ensure that diagonal element will be 1.
+                            if (j == i + 1)
+                            {
+                                mx = new Complex[,] { { 1 / a, 0.0 }, { 0.0, a } };
+                            }
+                        }
+                        else if (A[i, j - 1].Magnitude <= tol)
+                        {
+                            // Just swap columns with Pauli X matrix.
+                            mx = new Complex[,] { { 0.0, 1.0 }, { 1.0, 0.0 } };
+                            // But if it's last in a row, ensure that diagonal element will be 1.
+                            if (j == i + 1)
+                            {
+                                mx = new Complex[,] { { 0.0, b }, { 1 / b, 0.0 } };
+                            }
+                        }
+                        else
+                        {
+                            mx = MakeEliminatingMatrix(A[i, j - 1], A[i, j]);
+                        }
+                        var u2x2 = new TwoLevelUnitary(mx, j - 1, j);
+                        u2x2.MultiplyRight(A);
+                        u2x2.ConjugateTranspose();
+                        if (!u2x2.IsIdentity())
+                        {
+                            result.Add(u2x2);
+                        }
+                    }
+
+                    // At this point all elements in row i and in column i are zeros, with exception
+                    // of diagonal element A[i, i], which is equal to 1.
+                    Debug.Assert((A[i, i] - 1.0).Magnitude < tol);
+                }
+
+                var lastMatrix = new TwoLevelUnitary(new Complex[,] {
+                    { A[n - 2, n - 2], A[n - 2, n - 1] },
+                    { A[n - 1, n - 2], A[n - 1, n - 1] } }, n - 2, n - 1);
+                if (!lastMatrix.IsIdentity())
+                {
+                    result.Add(lastMatrix);
+                }
+                return result;
+            }
+
+            // Returns permutation of numbers from 0 to n-1 such that any two consequent numbers in
+            // it differ in exactly one bit.
+            private static int[] GrayCode(int n)
+            {
+                var result = new int[n];
+                for (int i = 0; i < n; i++)
+                {
+                    result[i] = i ^ (i / 2);
+                }
+                return result;
+            }
+
+            // Applies permutation to columns and rows of matrix.
+            private static Complex[,] PermuteMatrix(Complex[,] matrix, int[] perm)
+            {
+                int n = perm.Length;
+                Debug.Assert(matrix.GetLength(0) == n);
+                Debug.Assert(matrix.GetLength(1) == n);
+                var result = new Complex[n, n];
+                for (int i = 0; i < n; i++)
+                {
+                    for (int j = 0; j < n; j++)
+                    {
+                        result[i, j] = matrix[perm[i], perm[j]];
+                    }
+                }
+                return result;
+            }
+
+            // Returns list of two-level unitary matrices, which multiply to A.
+            // Every matrix has indices differing in exactly 1 bit.
+            private static IEnumerable<TwoLevelUnitary> TwoLevelDecomposeGray(Complex[,] A)
+            {
+                int n = A.GetLength(0);
+                Debug.Assert(A.GetLength(1) == n, "Matrix is not square.");
+                Debug.Assert(MatrixUtils.IsMatrixUnitary(A), "Matrix is not unitary.");
+                int[] perm = GrayCode(n);
+                A = PermuteMatrix(A, perm);
+                foreach (TwoLevelUnitary matrix in TwoLevelDecompose(A))
+                {
+                    matrix.ApplyPermutation(perm);
+                    matrix.OrderIndices();
+                    yield return matrix;
+                }
+            }
+
+            // Decomposes unitary matrix into product of 2-level unitary matrices.
+            // Every resulting 2-level matrix is represnted by 2x2 non-trivial submatrix and indices
+            // of this submatrix. Indices are guaranteed to be sorted and to differ in exactly one 
+            // bit.
+            private IQArray<(IQArray<IQArray<Quantum.Math.Complex>>, long, long)> Decompose(
+                IQArray<IQArray<Quantum.Math.Complex>> unitary)
+            {
+                Complex[,] a = MatrixUtils.MatrixFromQs(unitary);
+                return new QArray<(IQArray<IQArray<Quantum.Math.Complex>>, long, long)>(
+                    TwoLevelDecomposeGray(a).Select(matrix => matrix.ToQsharp())
+                );
+            }
+
+            public override Func<IQArray<IQArray<Quantum.Math.Complex>>,
+                IQArray<(IQArray<IQArray<Quantum.Math.Complex>>, long, long)>> __Body__ =>
+                Decompose;
+        }
+    }
+}