diff --git a/pyproject.toml b/pyproject.toml
index 27c379e..f47c40a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,7 +10,7 @@ name = "qmprs"
 version = "0.0.0rc0"
 dependencies = [
     "quick-core @ git+https://github.com/Qualition/quick",
-    "quimb == 1.10.0"
+    "quimb == 1.10.0",
 ]
 requires-python = ">=3.10, <3.13"
 authors = [
diff --git a/qmprs/primitives/mps.py b/qmprs/primitives/mps.py
index 0c75d9c..82b9c0b 100644
--- a/qmprs/primitives/mps.py
+++ b/qmprs/primitives/mps.py
@@ -48,45 +48,56 @@
 
 
 class MPS:
-    r""" `qmprs.primitives.MPS` is the class for creating and manipulating matrix product
-    states (MPS). This class wraps the `quimb.tensor.MatrixProductState` class to provide
-    a more user-friendly interface for creating and manipulating MPS.
+    r""" `qmprs.primitives.MPS` is the class for creating and manipulating matrix
+    product states (MPS). This class wraps the `quimb.tensor.MatrixProductState`
+    class to provide a more user-friendly interface for creating and manipulating
+    MPS.
 
     Refer to the link below for more information on the `quimb.tensor.MatrixProductState` class.
     https://quimb.readthedocs.io/en/latest/autoapi/quimb/tensor/tensor_1d/index.html#quimb.tensor.tensor_1d.MatrixProductState
 
     Notes
     -----
-    Matrix product states (MPS) are a class of 1D tensor networks that are widely used
-    in quantum computing to approximate the state of a quantum system. The MPS representation
-    allows for a polynomial or even exponential reduction in the number of parameters required
-    to represent a quantum state, making it a powerful tool for quantum state synthesis and
-    simulation.
-
-    The MPS representation is defined by first performing successive SVDs on the statevector
-    of the quantum system. We use SVD to find low-rank structure in the tensor network and
-    reduce the dimension of the tensors. We then choose a canonical form for the MPS, which
-    is either "left" or "right". The canonical form of the MPS states how the singular values
-    from the SVD are absorbed (contracted) with the left or right tensors.
-
-    The MPS can be further compressed by truncating the bond dimension of the MPS. The bond
-    dimension of the MPS determines the dimension of the unitary layers, and affects the
-    fidelity of the approximation. This is because the bond dimension captures the entanglement
-    structure of the quantum many-body system, where a higher bond dimension implies a higher
-    degree of entanglement.
-
-    Given each site will be represented as a $\chi\times\chi$ unitary matrix, the overall MPS
-    will have a scaling of $O(N\chi^2)$, where N is the number of sites and $\chi$ is the bond
-    dimension. Given a bond dimension of $2^{N/2}$ we can exactly represent any quantum state
-    of N qubits. However, for practical purposes, if we can keep the bond dimension constant,
-    the MPS will have a linear scaling with the number of sites.
+    Matrix product states (MPS) are a class of 1D tensor networks that are widely
+    used in quantum computing to approximate the state of a quantum system. The MPS
+    representation allows for a polynomial or even exponential reduction in the number
+    of parameters required to represent a quantum state, making it a powerful tool for
+    quantum state synthesis and simulation.
+
+    The MPS representation is defined by first performing successive SVDs on the
+    statevector of the quantum system. We use SVD to find low-rank structure in
+    the tensor network and reduce the dimension of the tensors. We then choose a
+    canonical form for the MPS, which is either "left" or "right". The canonical
+    form of the MPS states how the singular values from the SVD are absorbed
+    (contracted) with the left or right tensors.
+
+    The MPS can be further compressed by truncating the bond dimension of the MPS.
+    The bond dimension of the MPS determines the dimension of the unitary layers,
+    and affects the fidelity of the approximation. This is because the bond dimension
+    captures the entanglement structure of the quantum many-body system, where a higher
+    bond dimension implies a higher degree of entanglement.
+
+    MPS are particularly intended for approximating area-law entangled states, where
+    the entanglement entropy scales with the boundary of the system rather than the
+    volume. This is in contrast to states with volume-law entanglement, which cannot
+    be efficiently represented by MPS. This is a limitation of TNs in general, as
+    they require exponential bond dimension to represent such states. This is why
+    we opt for quantum computers to operate on such states.
+
+    Given each site will be represented as a $\chi\times\chi$ unitary matrix, the
+    overall MPS will have a scaling of $O(N\chi^2)$, where N is the number of sites
+    and $\chi$ is the bond dimension. Given a bond dimension of $2^{N/2}$ we can
+    exactly represent any quantum state of N qubits. However, for practical purposes,
+    if we can keep the bond dimension constant, the MPS will have a linear scaling
+    with the number of sites.
 
     The MPS can be written in the following
 
     $\ket{\psi} = \sum_{i_1, i_2, \cdots, i_N} Tr(A^{i_1}A^{i_2}\cdots A^{i_N}) \ket{i_1,i_2,\cdots,i_N}$
 
-    Where for arbitrary states, the MPS would be open-boundary condition, and non-translational
-    invariant, where A^{i} are not necessarily equal, and the first and last tensors are vectors.
+    Where for arbitrary states, the MPS would be open-boundary condition, and non-
+    translational invariant, where A^{i} are not necessarily equal, and the first
+    and last tensors are vectors.
 
     MPS Diagram:
     ```
@@ -96,11 +107,13 @@ class MPS:
     d d d d     d d d
     ```
 
-    where O represents the tensor at each site, d is the physical dimension, and D is the bond
-    dimension (also known as rank). For qubit systems, the physical dimension is 2.
+    where O represents the tensor at each site, d is the physical dimension, and
+    D is the bond dimension (also known as rank). For qubit systems, the physical
+    dimension is 2.
 
-    An important note is that the MPS representation is aimed for at least 2 qubits, as the MPS
-    approximates the entanglement structure of the quantum many-body systems.
+    An important note is that the MPS representation is aimed for at least 2 qubits,
+    as the MPS approximates the entanglement structure of the quantum many-body
+    systems.
 
     Parameters
     ----------
@@ -152,7 +165,7 @@ def __init__(
             self,
             statevector: Ket | NDArray[np.complex128] | None = None,
             mps: qtn.MatrixProductState | None = None,
-            bond_dimension: int=64
+            bond_dimension: int = 64
         ) -> None:
         """ Initialize a `qmprs.primitives.MPS` instance.
 
@@ -846,8 +859,16 @@ def generate_unitary_layer(self) -> UnitaryLayer:
 
         return generated_unitary_layer
 
-    def generate_bond_D_unitary_layer(self) -> UnitaryLayer:
-        r""" Truncate the unitary layer's bond dimension to 2.
+    def generate_bond_D_unitary_layer(
+            self,
+            optimize_truncated_mps: bool = False,
+            num_iterations_per_site: int = 25
+        ) -> UnitaryLayer:
+        r""" Truncate the unitary layer's bond dimension to 2. This method provides two
+        options: either to analytically compress the MPS to a bond dimension of 2, or to
+        variationally perform the fitting to make the compressed MPS as close as possible
+        to the original MPS. The latter is done by optimizing the compressed MPS to minimize
+        the infidelity with the original MPS.
 
         Notes
         -----
@@ -859,6 +880,14 @@ def generate_bond_D_unitary_layer(self) -> UnitaryLayer:
 
         https://arxiv.org/pdf/2209.00595, Figure 1
 
+        Parameters
+        ----------
+        `optimize_truncated_mps` : bool, optional, default=False
+            Whether to optimize the compressed MPS. If True, the compressed MPS will
+            be variationally optimized to minimize the infidelity with the original MPS.
+        `num_iterations_per_site` : int, optional, default=25
+            The number of iterations per site to optimize the compressed MPS.
+
         Returns
         -------
         `generated_unitary_layer` : UnitaryLayer
@@ -875,18 +904,39 @@ def generate_bond_D_unitary_layer(self) -> UnitaryLayer:
         """
         # Copy the MPS (as the MPS will be modified in place with
         # `.compress` and `.canonicalize` methods)
-        mps_copy = copy.deepcopy(self)
+        mps_truncated = copy.deepcopy(self)
 
         # Truncate the MPS to the bond dimension of 2 via SVD
-        mps_copy.compress(mode="right", max_bond_dimension=self.physical_dimension)
+        if optimize_truncated_mps:
+            max_iterations = num_iterations_per_site * self.num_sites
+
+            # using the 1-site method can be more efficient for fixed chi
+            # which we do via `bsz=1`
+            # The `permute_arrays` argument is used to ensure the
+            # resulting MPS is in the "lpr" form, which is what is used
+            # in `Sequential` encoding
+            # If you use different permute shape, and observe near 0
+            # fidelity, check to make sure this is set to the correct value
+            mps_truncated.mps = qtn.tensor_network_1d_compress(
+                mps_truncated.mps,
+                max_bond=2,
+                cutoff=0.0,
+                method="fit",
+                bsz=1,
+                max_iterations=max_iterations,
+                permute_arrays="lpr" # type: ignore
+            )
+            mps_truncated.bond_dimension = 2
+        else:
+            mps_truncated.compress(mode="right", max_bond_dimension=2)
 
         # To facilitate the loss-less conversion of all core
         # tensors in a TN into isometries (i.e. inner product
         # preserving transformations between Hilbert space)
         # we will canonicalize the MPS
-        mps_copy.canonicalize(mode="right", normalize=True)
+        mps_truncated.canonicalize(mode="right", normalize=True)
 
-        generated_unitary_layer = mps_copy.generate_unitary_layer()
+        generated_unitary_layer = mps_truncated.generate_unitary_layer()
 
         return generated_unitary_layer
 
@@ -973,7 +1023,7 @@ def _apply_inverse_unitary_layer(
     def apply_unitary_layer(
             self,
             unitary_layer: UnitaryLayer,
-            inverse: bool=False
+            inverse: bool = False
         ) -> None:
         """ Apply the unitary layer on the MPS. If inverse is True,
         we apply the inverse of the unitary layer to the MPS.
@@ -997,7 +1047,7 @@ def apply_unitary_layer(
     def apply_unitary_layers(
             self,
             unitary_layers: list[UnitaryLayer],
-            inverse: bool=False
+            inverse: bool = False
         ) -> None:
         """ Apply the unitary layers on the MPS. If inverse is True,
         we apply the inverse of the unitary layers in reverse order
@@ -1030,13 +1080,8 @@ def fidelity_with_zero_state(self) -> complex:
         -----
         >>> mps.fidelity_with_zero_state()
         """
-        zero_state = np.zeros(2**self.num_sites, dtype=np.complex128)
-        zero_state[0] = 1
-
-        # Compute the current statevector of the MPS
-        current_statevector = self.to_statevector(self.mps).data.flatten()
-
-        return np.dot(current_statevector.conj().T, zero_state)
+        zero_mps = qtn.MPS_computational_state([0] * self.num_sites, dtype='complex128')
+        return zero_mps @ self.mps # type: ignore
 
     def draw(self) -> plt.Figure:
         """ Draw the MPS.
@@ -1133,6 +1178,6 @@ def __eq__(
             value.mps.geometry_hash(strict_index_order=True)
 
         # Check if all the tensors in the MPSs are equal
-        all_close_eq = all(do("allclose", x, y) for x, y in zip(self.mps.arrays, self.mps.arrays))
+        all_close_eq = all(do("allclose", x, y) for x, y in zip(self.mps.arrays, value.mps.arrays))
 
         return geometry_hash_eq and all_close_eq
\ No newline at end of file
diff --git a/qmprs/synthesis/mps_encoding/base.py b/qmprs/synthesis/mps_encoding/base.py
index 3af2449..88ebb82 100644
--- a/qmprs/synthesis/mps_encoding/base.py
+++ b/qmprs/synthesis/mps_encoding/base.py
@@ -22,7 +22,7 @@
 from abc import ABC, abstractmethod
 import numpy as np
 from numpy.typing import NDArray
-from typing import Type, Literal
+from typing import Literal
 from quick.circuit import Circuit
 from quick.primitives import Ket
 
@@ -49,7 +49,7 @@ class MPSEncoder(ABC):
     """
     def __init__(
             self,
-            circuit_framework: Type[Circuit]
+            circuit_framework: type[Circuit]
         ) -> None:
         """ Initialize a `qmprs.mps_encoding.MPSEncoder` instance.
         """
@@ -59,7 +59,7 @@ def prepare_state(
             self,
             statevector: Ket | NDArray[np.complex128],
             bond_dimension: int,
-            compression_percentage: float=0.0,
+            compression_percentage: float = 0.0,
             index_type: Literal["row", "snake"]="row",
             **kwargs
         ) -> Circuit:
diff --git a/qmprs/synthesis/mps_encoding/sequential.py b/qmprs/synthesis/mps_encoding/sequential.py
index 34e58c5..595634f 100644
--- a/qmprs/synthesis/mps_encoding/sequential.py
+++ b/qmprs/synthesis/mps_encoding/sequential.py
@@ -24,7 +24,7 @@
 import quimb.tensor as qtn # type: ignore
 from quick.circuit import Circuit
 
-from qmprs.primitives.mps import MPS, UnitaryBlock, UnitaryLayer
+from qmprs.primitives.mps import MPS, UnitaryLayer
 from qmprs.synthesis.mps_encoding import MPSEncoder
 
 
@@ -35,16 +35,16 @@ class Sequential(MPSEncoder):
 
     We find the sequence of $\chi$ by $\chi$ unitary matrices that optimally
     disentangle the MPS to the product state $\ket{00\cdots 0}$. We compress
-    the max bond dimension to 2, so that we would only use one and two qubit gates.
-    We then reverse the sequence to obtain the quantum circuit that prepares
-    the MPS from the product state.
+    the max bond dimension to 2, so that we would only use one and two qubit
+    gates and reach $O(N)$ depth scaling. We then reverse the sequence to
+    obtain the quantum circuit that prepares the MPS from the product state.
 
     Each layer disentangles the MPS further, and depending on how entangled
     an MPS is, or how large, the number of layers needed to sufficiently
     disentangle the MPS will differ. The pseudo-code for the algorithm is
     available in [2] for developers' reference in Algorithm 1.
 
-    [1] Ran, Shi-Ju.
+    [1] Ran.
     Encoding of Matrix Product States into Quantum Circuits of One- and Two-Qubit Gates (2020).
     https://arxiv.org/abs/1908.07958
 
@@ -52,50 +52,86 @@ class Sequential(MPSEncoder):
     Decomposition of Matrix Product States into Shallow Quantum Circuits (2022).
     https://arxiv.org/abs/2209.00595
 
+    [3] Lin, Dilip, Green, Smith, Pollmann.
+    Real- and imaginary-time evolution with compressed quantum circuits (2008).
+    https://arxiv.org/pdf/2008.10322
+
     Notes
     -----
     The sequential encoding is a method to prepare a target MPS using a sequence
     of $\chi x \chi$ unitary matrices on each site (aka qubit), where $\chi$ is
     the bond dimension of the MPS.
 
-    Assuming the bond dimension is a power of two, we would require $n = \log_2(\chi)$
-    qubits to prepare the unitary matrix. The exact preparation of a general unitary
-    matrix scales $O(2^{2n})$, and it is equivalent to $O(\chi^2)$. Furthermore,
-    given the linear dependence of the unitary matrices on the number of sites,
-    the overall scaling is $O(N\chi^2)$ as stated in [1].
-
-    In the implementation of the sequential encoding, we provide two parameters:
-    1) `num_layers` : The number of unitary layers used to prepare the MPS.
-    2) `bond_dimension` : The maximum bond dimension of the MPS.
-
-    Unlike [1], the bond dimension does not directly influence the circuit depth.
-    Instead, the bond dimension controls the maximum possible fidelity achievable
-    through the MPS encoding. The number of layers is the primary parameter that
-    influences the circuit depth.
-
-    The algorithm was designed for long-range correlation and follows a unitary-only
-    operation approach which limits the efficiency of the encoding. The number
-    of layers scale linearly with the circuit depth, and exponentially with the
-    number of sites. However, the exponential growth is significantly reduced by
-    the low-rank structure of the MPS representation, which allows for increasingly
-    efficient encoding of quantum states as we scale the number of sites when compared
-    to exact encoding schema such as Mottonen, Shende, or SOTA Isometry by Iten et al.
-    Additionally, given the analytical decomposition employed, the sequential encoding
-    is computationally more efficient, however, that also means that increasing the
-    number of layers will only slightly improve the fidelity of the encoding.
+    Assuming the bond dimension is a power of two, we would require
+    $n = \log_2(\chi) + 1$ qubits to prepare each tensor depending on its bond
+    dimension. The exact preparation of a general unitary matrix scales $O(2^{2n})$,
+    and it is equivalent to $O(\chi^2)$.
+
+    Furthermore, given the linear dependence of the unitary matrices on the number
+    of sites, the overall scaling is $O(N\chi^2)$. We can exactly prepare the MPS
+    using a single layer but due to the impact of bond dimension on the unitary sizes,
+    we would require multi-qubit gates for certain tensors ([3] Appendix A).
+
+    What Ran [1] proposes instead is to truncate the bond dimension of the MPS to 2,
+    and prepare the truncated MPS instead. This keeps the bond dimension constant,
+    thus $O(N)$ scaling. Given truncation causes loss of fidelity, we need to use
+    multiple layers instead to approximately prepare the target MPS.
+
+    Given the algorithm utilizes bond 2 truncation, it is intended for area-law
+    entangled states which do not require exponential bond dimensions to represent,
+    and thus do not lose significant fidelity when truncated to bond 2. This also
+    means significantly fewer layers needed to prepare. The synthesis may work for
+    volume-law entangled states, but it is not guaranteed to be efficient or
+    effective.
+
+    Ran's approach analytically produces the unitary layers that disentangle the MPS
+    and is known to slowly converge to the product state. This makes it increasingly
+    harder to reach a target fidelity as the fidelity gain slows down and plateaus
+    with each layers.
 
     To achieve a higher fidelity within a reasonable circuit depth, we use environment
-    tensor updates to reach the optimal gates for the circuit based on [2].
+    tensor updates to reach the optimal gates for the circuit based on [2]. This approach
+    allows us to efficiently optimize the unitary layers by sweeping through the tensor
+    network representation of the circuit (produced from the unitary layers) and updating
+    the tensors to maximize the inner product with the target MPS.
+
+    Additionally, to further improve the fidelity for volume-law entangled states, we
+    variationally optimize the bond 2 truncation of the MPS to improve the fidelity
+    between the truncated MPS and the target MPS. This provides slight improvement in
+    the fidelity for area-law entangled states, but is considerably more effective for
+    volume-law entangled states when paired with environment tensor updates. Due to the
+    stochastic nature of the optimization, one may need to re-run the compilation multiple
+    times to get the best fidelity.
 
     Parameters
     ----------
     `circuit_framework` : type[quick.circuit.Circuit]
         The quantum circuit framework.
+    `variationally_optimize_truncated_mps` : bool, optional, default=False
+        Whether to variationally optimize the bond 2 truncation of the MPS to improve
+        the fidelity between the truncated MPS and the target MPS. This is useful for
+        volume-law entangled states where the bond dimension is significantly larger
+        than 2, and the truncation causes significant loss of fidelity. This is turned
+        off by default, as it slightly increases the runtime of the compilation and
+        provides only marginal improvement in fidelity for area-law entangled states.
+        Furthermore, given the stochastic nature of the optimization, one may need to
+        re-run the compilation multiple times to get the best fidelity, and thus it may
+        affect the reproducibility and in turn cause some unit tests to fail.
+    `num_iterations_per_site` : int, optional, default=25
+        The number of iterations to perform for the variational optimization of the
+        bond 2 truncation of the MPS. This defines `max_iterations` for the quimb
+        `qtn.tensor_network_1d_compress` method via `num_iterations_per_site * mps.num_sites`.
 
     Attributes
     ----------
     `circuit_framework` : type[quick.circuit.Circuit]
         The quantum circuit framework.
+    `variationally_optimize_truncated_mps` : bool
+        Whether to variationally optimize the bond 2 truncation of the MPS to improve
+        the fidelity between the truncated MPS and the target MPS.
+    `num_iterations_per_site` : int
+        The number of iterations to perform for the variational optimization of the
+        bond 2 truncation of the MPS.
     `fidelity_threshold` : float, optional, default=0.999999
         The fidelity threshold for the MPS encoding. The encoding stops when the
         fidelity of the MPS with the product state is greater than or equal to the
@@ -104,7 +140,8 @@ class Sequential(MPSEncoder):
     Raises
     ------
     ValueError
-        If the number of layers is not a positive integer.
+        - If the number of layers is not a positive integer.
+        - If the number of sweeps is not a non-negative integer.
 
     Usage
     -----
@@ -112,11 +149,16 @@ class Sequential(MPSEncoder):
     """
     def __init__(
             self,
-            circuit_framework: type[Circuit]
+            circuit_framework: type[Circuit],
+            variationally_optimize_truncated_mps: bool = False,
+            num_iterations_per_site: int = 25
         ) -> None:
 
         super().__init__(circuit_framework)
 
+        self.variationally_optimize_truncated_mps = variationally_optimize_truncated_mps
+        self.num_iterations_per_site = num_iterations_per_site
+
         self._fidelity_threshold = 1 - 1e-6
 
     @property
@@ -155,7 +197,7 @@ def fidelity_threshold(
     @staticmethod
     def _apply_unitary_layer_to_circuit(
             circuit: Circuit,
-            unitary_layer: list[UnitaryBlock]
+            unitary_layer: UnitaryLayer
         ) -> None:
         """ Apply a unitary layer to the quantum circuit.
 
@@ -163,7 +205,7 @@ def _apply_unitary_layer_to_circuit(
         ----------
         `circuit` : quick.circuit.Circuit
             The quantum circuit.
-        `unitary_layer` : list[qtn.Tensor]
+        `unitary_layer` : UnitaryLayer
             The unitary layer to be applied to the circuit.
         """
         for start_index, end_index, unitary_block in unitary_layer:
@@ -317,7 +359,10 @@ def _get_unitary_layer(
         """
         # Generate the bond 2 truncation of the unitary layer
         # to form one and two qubit gates given Fig. 1 in [1]
-        unitary_layer = mps.generate_bond_D_unitary_layer()
+        unitary_layer = mps.generate_bond_D_unitary_layer(
+            optimize_truncated_mps=self.variationally_optimize_truncated_mps,
+            num_iterations_per_site=self.num_iterations_per_site
+        )
 
         # Given MPS ~= U_k|00...0>, we need to apply the inverse of U_k
         # to disentangle the MPS to the product state |00...0>
@@ -407,7 +452,7 @@ def _sweep_unitary_layers(
         using the environment tensor updates to reach the optimal gates
         for the circuit based on [2].
 
-        This method implements Oall and allows for Iter Oi from [2].
+        Synonymously, this method implements Oall from [2].
 
         Notes
         -----
@@ -596,5 +641,7 @@ def prepare_mps(
 
         if not isinstance(num_layers, int) or num_layers < 1:
             raise ValueError("The number of layers must be a positive integer.")
+        if not isinstance(num_sweesps, int) or num_sweesps < 0:
+            raise ValueError("The number of sweeps must be a non-negative integer.")
 
         return self._sequential_unitary_circuit(mps, num_layers, num_sweesps)
\ No newline at end of file
diff --git a/tests/synthesis/mps_encoding/test_sequential_encoding.py b/tests/synthesis/mps_encoding/test_sequential_encoding.py
index 22f103c..299c691 100644
--- a/tests/synthesis/mps_encoding/test_sequential_encoding.py
+++ b/tests/synthesis/mps_encoding/test_sequential_encoding.py
@@ -18,6 +18,7 @@
 
 import numpy as np
 from numpy.typing import NDArray
+import pytest
 from quick.circuit import QiskitCircuit
 
 from qmprs.primitives import MPS
@@ -43,6 +44,32 @@ def generate_random_state(num_qubits: int) -> NDArray[np.complex128]:
     statevector /= np.linalg.norm(statevector)
     return statevector
 
+def generate_random_clifford_circuit(num_qubits: int) -> NDArray[np.complex128]:
+    """ Generate a random Clifford circuit state.
+
+    Parameters
+    ----------
+    num_qubits : int
+        The number of qubits.
+
+    Returns
+    -------
+    NDArray[np.complex128]
+        The random Clifford circuit state.
+    """
+    from qiskit.circuit.random import random_clifford_circuit
+    from qiskit.quantum_info import Statevector
+
+    gates = ["cx", "cz", "cy", "swap", "x", "y", "z", "s", "sdg", "h"]
+    qc = random_clifford_circuit(
+        num_qubits,
+        gates=gates, # type: ignore
+        num_gates=10 * num_qubits * num_qubits,
+        seed=1,
+    )
+
+    return Statevector(qc).data
+
 
 class TestSequential(Template):
     """ `tests.synthesis.mps_encoding.TestSequential` is the tester for `qmprs.synthesis.mps_encoding.Sequential` class.
@@ -58,13 +85,17 @@ def test_prepare_state(self) -> None:
         encoder = Sequential(circuit_framework=QiskitCircuit)
 
         # Prepare the MPS from the statevector using the Sequential encoder
-        circuit = encoder.prepare_state(statevector=statevector, bond_dimension=32, num_layers=32)
+        circuit = encoder.prepare_state(
+            statevector=statevector,
+            bond_dimension=32,
+            num_layers=32
+        )
 
         # Extract the statevector from the circuit
         statevector_from_circuit = circuit.get_statevector()
 
         # Ensure that the statevector from the circuit is equal to the original statevector
-        assert 1 - abs(np.dot(statevector_from_circuit.conj(), statevector)) < 1e-2
+        assert 1 - abs(np.vdot(statevector_from_circuit, statevector)) < 1e-2
 
     def test_prepare_mps(self) -> None:
         """ Test the preparation of the MPS from a MPS.
@@ -86,7 +117,7 @@ def test_prepare_mps(self) -> None:
         statevector_from_circuit = circuit.get_statevector()
 
         # Ensure that the statevector from the circuit is equal to the original statevector
-        assert 1 - abs(np.dot(statevector_from_circuit.conj(), statevector)) < 1e-2
+        assert 1 - abs(np.vdot(statevector_from_circuit, statevector)) < 1e-2
 
     def test_prepare_circuit_with_partial_entanglement(self) -> None:
         """ Test the preparation of the MPS from a statevector with partial entanglement.
@@ -107,13 +138,17 @@ def test_prepare_circuit_with_partial_entanglement(self) -> None:
         encoder = Sequential(circuit_framework=QiskitCircuit)
 
         # Prepare the MPS from the statevector using the Sequential encoder
-        circuit = encoder.prepare_state(statevector=statevector, bond_dimension=32, num_layers=1)
+        circuit = encoder.prepare_state(
+            statevector=statevector,
+            bond_dimension=32,
+            num_layers=1
+        )
 
         # Extract the statevector from the circuit
         statevector_from_circuit = circuit.get_statevector()
 
         # Ensure that the statevector from the circuit is equal to the original statevector
-        assert 1 - abs(np.dot(statevector_from_circuit.conj(), statevector)) < 1e-2
+        assert 1 - abs(np.vdot(statevector_from_circuit, statevector)) < 1e-2
 
         # The produced circuit is much shallower compared to a full entangled circuit
         # and we need to check that the circuit depth is less than 20
@@ -140,14 +175,17 @@ def test_prepare_circuit_with_partial_entanglement_with_sweep(self) -> None:
 
         # Prepare the MPS from the statevector using the Sequential encoder
         circuit = encoder.prepare_state(
-            statevector=statevector, bond_dimension=32, num_layers=1, num_sweeps=1
+            statevector=statevector,
+            bond_dimension=32,
+            num_layers=1,
+            num_sweeps=1
         )
 
         # Extract the statevector from the circuit
         statevector_from_circuit = circuit.get_statevector()
 
         # Ensure that the statevector from the circuit is equal to the original statevector
-        assert 1 - abs(np.dot(statevector_from_circuit.conj(), statevector)) < 1e-2
+        assert 1 - abs(np.vdot(statevector_from_circuit, statevector)) < 1e-2
 
         # The produced circuit is much shallower compared to a full entangled circuit
         # and we need to check that the circuit depth is less than 20
@@ -164,17 +202,22 @@ def test_layer_improvement(self) -> None:
         # Define the Sequential encoder
         encoder = Sequential(circuit_framework=QiskitCircuit)
 
-        layer_fidelity = []
+        layer_fidelity: list[float] = []
 
         for i in range(1, 10):
             # Prepare the MPS from the statevector using the Sequential encoder
-            circuit = encoder.prepare_state(statevector=statevector, bond_dimension=64, num_layers=i)
+            circuit = encoder.prepare_state(
+                statevector=statevector,
+                bond_dimension=64,
+                num_layers=i
+            )
 
             # Extract the statevector from the circuit
             statevector_from_circuit = circuit.get_statevector()
 
-            # Compute the fidelity between the statevector from the circuit and the original statevector
-            fidelity = abs(np.dot(statevector_from_circuit.conj(), statevector))
+            # Compute the fidelity between the statevector from the circuit and the
+            # original statevector
+            fidelity = float(abs(np.vdot(statevector_from_circuit, statevector)))
             layer_fidelity.append(fidelity)
 
         # Ensure that the fidelity increases with the number of layers
@@ -190,18 +233,86 @@ def test_sweep_improvement(self) -> None:
         # Define the Sequential encoder
         encoder = Sequential(circuit_framework=QiskitCircuit)
 
-        bond_fidelity = []
+        bond_fidelity: list[float] = []
 
         for i in range(1, 8):
             # Prepare the MPS from the statevector using the Sequential encoder
-            circuit = encoder.prepare_state(statevector=statevector, bond_dimension=64, num_layers=6, num_sweeps=i)
+            circuit = encoder.prepare_state(
+                statevector=statevector,
+                bond_dimension=64,
+                num_layers=6,
+                num_sweeps=i
+            )
 
             # Extract the statevector from the circuit
             statevector_from_circuit = circuit.get_statevector()
 
-            # Compute the fidelity between the statevector from the circuit and the original statevector
-            fidelity = abs(np.dot(statevector_from_circuit.conj(), statevector))
+            # Compute the fidelity between the statevector from the circuit and the
+            # original statevector
+            fidelity = float(abs(np.vdot(statevector_from_circuit, statevector)))
             bond_fidelity.append(fidelity)
 
         # Ensure that the fidelity increases with the bond dimension
-        assert np.all(np.diff(bond_fidelity) >= 0)
\ No newline at end of file
+        assert np.all(np.diff(bond_fidelity) >= 0)
+
+    @pytest.mark.parametrize("num_qubits", [8, 10, 12, 14])
+    def test_optimize_bond_2_truncation(
+            self,
+            num_qubits: int
+        ) -> None:
+        """ Test the optimization of the bond 2 truncation.
+
+        Parameters
+        ----------
+        `num_qubits` : int
+            The number of qubits.
+        """
+        # Define the number of qubits and generate a random statevector
+        statevector = generate_random_clifford_circuit(num_qubits)
+
+        # Define the Sequential encoder
+        encoder_without_optimization = Sequential(circuit_framework=QiskitCircuit)
+        encoder_with_optimization = Sequential(
+            circuit_framework=QiskitCircuit,
+            variationally_optimize_truncated_mps=True,
+            num_iterations_per_site=15
+        )
+
+        # Prepare the MPS from the statevector using the Sequential encoder
+        # without optimization
+        non_optimized_circuit = encoder_without_optimization.prepare_state(
+            statevector=statevector,
+            bond_dimension=2**num_qubits,
+            num_layers=3,
+            num_sweeps=50
+        )
+        non_optimized_statevector = non_optimized_circuit.get_statevector()
+        non_optimized_fidelity = np.vdot(statevector, non_optimized_statevector)
+
+        best_circuit = None
+        best_fidelity = 0.0
+
+        # Prepare the MPS from the statevector using the Sequential encoder
+        # with optimization
+        # Given that the optimization is stochastic, we will run it
+        # multiple times and keep the best result
+        for _ in range(10):
+            optimized_circuit = encoder_with_optimization.prepare_state(
+                statevector=statevector,
+                bond_dimension=2**num_qubits,
+                num_layers=3,
+                num_sweeps=50
+            )
+
+            fidelity = np.vdot(statevector, optimized_circuit.get_statevector())
+
+            if best_circuit is None:
+                best_circuit = optimized_circuit
+                best_fidelity = fidelity
+            elif fidelity > best_fidelity:
+                best_circuit = optimized_circuit
+                best_fidelity = fidelity
+
+        # Ensure that the statevector from the optimized circuit is closer to the
+        # original statevector than the non-optimized circuit
+        assert abs(best_fidelity) > abs(non_optimized_fidelity)
\ No newline at end of file