feature: MPS Encoding Algorithm

README.md

@@ -18,6 +18,7 @@ Running notebooks locally requires additional dependencies located in [notebooks
 | CHSH Inequality | [CHSH_Inequality.ipynb](notebooks/textbook/CHSH_Inequality.ipynb) | [Clauser1970](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.23.880) |
 | Deutsch-Jozsa | [Deutsch_Jozsa_Algorithm.ipynb](notebooks/textbook/Deutsch_Jozsa_Algorithm.ipynb) | [Deutsch1992](https://royalsocietypublishing.org/doi/10.1098/rspa.1992.0167) |
 | Grover's Search | [Grovers_Search.ipynb](notebooks/textbook/Grovers_Search.ipynb) | [Figgatt2017](https://www.nature.com/articles/s41467-017-01904-7), [Baker2019](https://arxiv.org/abs/1904.01671) |
+| MPS Encoding | [MPS_Encoding.ipynb](notebooks/textbook/MPS_Encoding.ipynb) | [Ran2019](https://arxiv.org/abs/1908.07958), [Rudolph2022](https://arxiv.org/abs/2209.00595) |
 | QAOA | [Quantum_Approximate_Optimization_Algorithm.ipynb](notebooks/textbook/Quantum_Approximate_Optimization_Algorithm.ipynb) | [Farhi2014](https://arxiv.org/abs/1411.4028) |
 | Quantum Circuit Born Machine | [Quantum_Circuit_Born_Machine.ipynb](notebooks/textbook/Quantum_Circuit_Born_Machine.ipynb) | [Benedetti2019](https://www.nature.com/articles/s41534-019-0157-8),  [Liu2018](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.98.062324) | 
 | QFT | [Quantum_Fourier_Transform.ipynb](notebooks/textbook/Quantum_Fourier_Transform.ipynb) | [Coppersmith2002](https://arxiv.org/abs/quant-ph/0201067) | 

setup.py

@@ -34,6 +34,8 @@
         "scipy>=1.5.2",
         # Sympy 1.13 produces different results for Simon's algorithm
         "sympy<1.13",
+        "quimb",
+        "qiskit_braket_provider",
     ],
     extras_require={
         "test": [

src/braket/experimental/algorithms/mps_encoding/__init__.py

@@ -0,0 +1,16 @@
+# Copyright Amazon.com Inc. or its affiliates. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License"). You
+# may not use this file except in compliance with the License. A copy of
+# the License is located at
+#
+#     http://aws.amazon.com/apache2.0/
+#
+# or in the "license" file accompanying this file. This file is
+# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF
+# ANY KIND, either express or implied. See the License for the specific
+# language governing permissions and limitations under the License.
+
+from braket.experimental.algorithms.mps_encoding.mps_encoding import (  # noqa: F401
+    mps_circuit_from_statevector,
+)

src/braket/experimental/algorithms/mps_encoding/mps_encoding.md

@@ -0,0 +1,10 @@
+Matrix Product State (MPS) encoding tackles the critical task of state preparation in $O(N)$ circuit depth, an exponential reduction in circuit depth compared to
+exact encoding algorithms like Shende, Isometry, and Mottonen. The current implementation uses an analytical decomposition based on work by
+Ran et al. to create a staircase ansatz comprised of 1 and 2 qubit unitary gates. Main consideration with the decomposition is that it is prone to
+approximation limitations imposed by MPS IR, which limits the adequate performance of the synthesis to area-law entangled states only.
+
+<!--
+[metadata-name]: MPS Encoding
+[metadata-tags]: Textbook
+[metadata-url]: https://github.com/amazon-braket/amazon-braket-algorithm-library/tree/main/src/braket/experimental/algorithms/mps_encoding
+-->

src/braket/experimental/algorithms/mps_encoding/mps_encoding.py

@@ -0,0 +1,316 @@
+# Copyright Amazon.com Inc. or its affiliates. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License"). You
+# may not use this file except in compliance with the License. A copy of
+# the License is located at
+#
+#     http://aws.amazon.com/apache2.0/
+#
+# or in the "license" file accompanying this file. This file is
+# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF
+# ANY KIND, either express or implied. See the License for the specific
+# language governing permissions and limitations under the License.
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import numpy as np
+import quimb.tensor as qtn
+from qiskit_braket_provider import to_braket
+
+from braket.circuits import Circuit
+
+if TYPE_CHECKING:
+    import numpy.typing as npt
+
+
+def _gram_schmidt(matrix: npt.NDArray[np.complex128]) -> npt.NDArray[np.complex128]:
+    """Perform Gram-Schmidt orthogonalization on the columns of a square matrix
+    to define the unitary block to encode the MPS.
+
+    This implementation provides slightly better fidelity compared to the Schmidt
+    decomposition approach used in scipy.linalg.qr for our use case.
+
+    Args:
+        matrix (NDArray[np.complex128]): Square matrix with complex entries.
+
+    Returns:
+        (NDArray[np.complex128]): A unitary matrix with orthonormal columns derived from the input matrix.
+        If a column is approximately zero, it is replaced with a random vector.
+    """
+    assert len(matrix.shape) == 2 and matrix.shape[0] == matrix.shape[1], (
+        "Input matrix must be square."
+    )
+
+    num_rows: int = matrix.shape[0]
+    unitary = np.zeros((num_rows, num_rows), dtype=np.complex128)
+    orthonormal_basis: list[npt.NDArray[np.complex128]] = []
+
+    for j in range(num_rows):
+        column = matrix[:, j]
+
+        # If column is (approximately) zero, replace with random
+        if np.allclose(column, 0):
+            column = np.random.uniform(-1, 1, num_rows) + 1j * np.random.uniform(-1, 1, num_rows)
+
+        # Gram-Schmidt orthogonalization
+        for basis_vector in orthonormal_basis:
+            column -= (basis_vector.conj().T @ column) * basis_vector
+
+        # Handle near-zero vectors (linear dependence)
+        if np.linalg.norm(column) < 1e-12:  # pragma: no cover
+            column = np.random.uniform(-1, 1, num_rows) + 1j * np.random.uniform(-1, 1, num_rows)
+
+            for basis_vector in orthonormal_basis:
+                column -= (basis_vector.conj().T @ column) * basis_vector
+
+        unitary[:, j] = column / np.linalg.norm(column)
+        orthonormal_basis.append(unitary[:, j])
+
+    return unitary
+
+
+class MPSSequential:
+    """MPS to Circuit synthesis via sequential layers of unitaries.
+
+    This class provides methods to convert a Matrix Product State (MPS)
+    into a quantum circuit representation using layers of two-qubit and single-qubit unitaries
+    in the form of a staircase.
+
+    This approach allows for approximately encoding an area-law entangled state
+    using O(N) circuit depth, where N is the number of qubits by leveraging the MPS structure
+    as an intermediate representation. This is an exponential improvement over the
+    exact method of encoding arbitrary states, which typically requires O(2^N) depth.
+
+    The analytical decomposition is based on the approach presented in:
+    [arXiv:1908.07958](https://arxiv.org/abs/1908.07958){:.external}
+    """
+
+    def __init__(self, max_fidelity_threshold: float = 0.95) -> None:
+        """Initialize the Sequential class.
+
+        Args:
+            max_fidelity_threshold (float): The maximum fidelity required, after
+                which we can stop the encoding to save depth.
+        """
+        self.max_fidelity_threshold = max_fidelity_threshold
+
+    def generate_layer(
+        self, mps: qtn.MatrixProductState
+    ) -> list[tuple[list[int], npt.NDArray[np.complex128]]]:
+        """Convert a Matrix Product State (MPS) to a circuit representation
+        using a single unitary layer.
+
+        Args:
+            mps (qtn.MatrixProductState): The MPS to convert.
+
+        Returns:
+            (list[tuple[list[int], npt.NDArray[np.complex128]]]): A list of tuples
+                representing the unitary layer of the circuit.
+                Each tuple contains:
+                - A list of qubit indices (in LSB order) that the unitary acts on.
+                - A unitary matrix (as a 2D NumPy array) that encodes the MPS.
+        """
+        num_sites = mps.L
+
+        unitary_layer: list[tuple[list[int], npt.NDArray[np.complex128]]] = []
+
+        for i, tensor in enumerate(reversed(mps.arrays)):
+            i = num_sites - i - 1
+
+            # MPS representation uses 1D entanglement, thus we need to define
+            # the range of the indices via the tensor shape
+            # i.e., if q0 and q3 are entangled, then regardless of q1 and q2 being
+            # entangled the entanglement range would be q0-q3
+            if i == 0:
+                d_right, d = tensor.shape
+                tensor = tensor.reshape((1, d_right, d))
+            if i == num_sites - 1:
+                d_left, d = tensor.shape
+                tensor = tensor.reshape((d_left, 1, d))
+
+            tensor = np.swapaxes(tensor, 1, 2)
+
+            # Combine the physical index and right-virtual index of the tensor to construct
+            # an isometry matrix
+            d_left, d, d_right = tensor.shape
+            isometry = tensor.reshape((d * d_left, d_right))
+
+            qubits = [
+                num_sites - 1 - q for q in range(i, i - int(np.ceil(np.log2(d_left))) - 1, -1)
+            ]
+
+            matrix = np.zeros((isometry.shape[0], isometry.shape[0]), dtype=isometry.dtype)
+
+            # Keep columns for which all ancillas are in the zero state
+            matrix[:, : isometry.shape[1]] = isometry
+
+            # Perform Gram-Schmidt orthogonalization to ensure the columns are orthonormal
+            unitary = _gram_schmidt(matrix)
+
+            unitary_layer.append((qubits, unitary))
+
+        return unitary_layer
+
+    def mps_to_circuit_approx(
+        self, statevector: npt.NDArray[np.complex128], max_num_layers: int, chi_max: int
+    ) -> Circuit:
+        """Approximately encodes the MPS into a circuit via multiple layers
+        of exact encoding of bond 2 truncated MPS.
+
+        Whilst we can encode the MPS exactly in a single layer, we require
+        $log(chi) + 1$ qubits for each tensor, which results in larger circuits.
+        This function uses bond 2 which allows us to encode the MPS using one and
+        two qubit gates, which results in smaller circuits, and easier to run on
+        hardware.
+
+        This is the core idea of Ran's paper [1].
+
+        [1] https://arxiv.org/abs/1908.07958
+
+        Args:
+            statevector (NDArray[np.complex128]): The statevector to convert.
+            max_num_layers (int): The number of layers to use in the circuit.
+            chi_max (int): The maximum bond dimension of the target MPS.
+
+        Returns:
+            (Circuit): The generated quantum circuit that encodes the MPS.
+        """
+        mps_dense = qtn.MatrixProductState.from_dense(statevector)
+        mps: qtn.MatrixProductState = qtn.tensor_1d_compress.tensor_network_1d_compress(
+            mps_dense, max_bond=chi_max
+        )
+        mps.permute_arrays()
+
+        mps.compress(form="left", max_bond=chi_max)
+        mps.left_canonicalize(normalize=True)
+
+        compressed_mps = mps.copy(deep=True)
+        disentangled_mps = mps.copy(deep=True)
+
+        circuit = Circuit()
+
+        unitary_layers = []
+
+        # Initialize the zero state |00...0> to serve as comparison
+        # for the disentangled MPS
+        zero_state = np.zeros((2**mps.L,), dtype=np.complex128)
+        zero_state[0] = 1.0
+
+        # Ran's approach uses a iterative disentanglement of the MPS
+        # where each layer compresses the MPS to a maximum bond dimension of 2
+        # and applies the inverse of the layer to disentangle the MPS
+        # After a few layers we are adequately close to |00...0> state
+        # after which we can simply reverse the layers (no inverse) and apply them
+        # to the |00...0> state to obtain the MPS state
+        for _ in range(max_num_layers):
+            # Compress the MPS from the previous layer to a maximum bond dimension of 2,
+            # |ψ_k> -> |ψ'_k>
+            compressed_mps = disentangled_mps.copy(deep=True)
+
+            # Normalization improves fidelity of the encoding
+            compressed_mps.normalize()
+            compressed_mps.compress(form="left", max_bond=2)
+
+            unitary_layers.append(self.generate_layer(compressed_mps))
+
+            # To update the MPS definition, apply the inverse of U_k to disentangle |ψ_k>,
+            # |ψ_(k+1)> = inv(U_k) @ |ψ_k>
+            for i, (_, U) in enumerate(reversed(unitary_layers[-1])):
+                inverse = U.conj().T
+
+                if inverse.shape[0] == 4:
+                    disentangled_mps.gate_split_(inverse, (i - 1, i))
+                else:
+                    disentangled_mps.gate_(inverse, (i), contract=True)
+
+            # Compress the disentangled MPS to a maximum bond dimension of chi_max
+            # This is to ensure that the disentangled MPS does not grow too large
+            # and improves the fidelity of the encoding
+            disentangled_mps: qtn.MatrixProductState = (  # type: ignore
+                qtn.tensor_1d_compress.tensor_network_1d_compress(
+                    disentangled_mps, max_bond=chi_max
+                )
+            )
+
+            fidelity = np.abs(np.vdot(disentangled_mps.to_dense(), zero_state)) ** 2
+            fidelity = float(np.real(fidelity))
+
+            if fidelity >= self.max_fidelity_threshold:
+                break
+
+        # The layers disentangle the MPS to a state close to |00...0>
+        # inv(U_k) ... inv(U_1) |ψ> = |00...0>
+        # So, we have to reverse the layers and apply them to the |00...0> state
+        # to obtain the MPS state
+        # |ψ> = U_1 ... U_k |00...0>
+        unitary_layers.reverse()
+
+        for unitary_layer in unitary_layers:
+            for qubits, unitary in unitary_layer:
+                qubits = list(reversed([mps.L - 1 - q for q in qubits]))
+
+                # In certain cases, floating point errors can cause `unitary`
+                # to not be unitary
+                # We handle this by using SVD to approximate it again with the
+                # closest unitary
+                U, _, Vh = np.linalg.svd(unitary)
+                unitary = U @ Vh
+
+                circuit.unitary(matrix=unitary, targets=qubits)
+
+        return to_braket(circuit, add_measurements=False, basis_gates=["u3", "cx", "global_phase"])
+
+    def __call__(
+        self, statevector: npt.NDArray[np.complex128], max_num_layers: int = 10
+    ) -> Circuit:
+        """Call the instance to create the circuit that encodes the statevector.
+
+        Args:
+            statevector (NDArray[np.complex128]): The statevector to convert.
+
+        Returns:
+            (Circuit): The generated quantum circuit.
+        """
+        num_qubits = int(np.ceil(np.log2(len(statevector))))
+
+        # Single qubit statevector is optimal, and cannot be
+        # further improved given depth of 1
+        if num_qubits == 1:
+            circuit = Circuit()
+
+            magnitude = abs(statevector)
+            phase = np.angle(statevector)
+
+            divisor = magnitude[0] ** 2 + magnitude[1] ** 2
+            alpha_y = 2 * np.arcsin(np.sqrt(magnitude[1] ** 2 / divisor)) if divisor != 0 else 0.0
+            alpha_z = phase[1] - phase[0]
+            global_phase = sum(phase / len(statevector))
+
+            circuit.ry(angle=alpha_y, target=0)
+            circuit.rz(angle=alpha_z, target=0)
+            circuit.gphase(angle=global_phase)
+
+            return circuit
+
+        circuit = self.mps_to_circuit_approx(statevector, max_num_layers, 2**num_qubits)
+
+        return circuit
+
+
+def mps_circuit_from_statevector(
+    statevector: npt.NDArray[np.complex128], max_num_layers: int = 10, target_fidelity: float = 0.95
+) -> Circuit:
+    """Create the circuit that encodes the statevector using MPS synthesis.
+
+    Args:
+        statevector (NDArray[np.complex128]): The target statevector to be encoded.
+        max_num_layers (int): The maximum number of layers allowed in the circuit.
+        target_fidelity (float): The target fidelity for the approximation.
+
+    Returns:
+        A braket.Circuit that encodes the statevector.
+    """
+    encoder = MPSSequential(max_fidelity_threshold=target_fidelity)
+    return encoder(statevector, max_num_layers=max_num_layers)