diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 068d357c9c0..cb5f4a98005 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -186,10 +186,15 @@
   interface-specific scalar data, eg `[(tf.Variable(1.1), tf.Variable(2.2))]`.
   [(#4603)](https://github.com/PennyLaneAI/pennylane/pull/4603)
 
+* When decomposing a unitary matrix with `one_qubit_decomposition`, and opting to include the `GlobalPhase` 
+  in the decomposition, the phase is no longer cast to `dtype=complex`.
+  [(#4653)](https://github.com/PennyLaneAI/pennylane/pull/4653)
+
 * `qml.cut_circuit` is now compatible with circuits that compute the expectation values of Hamiltonians 
   with two or more terms.
   [(#4642)](https://github.com/PennyLaneAI/pennylane/pull/4642)
 
+
 * `_qfunc_output` has been removed from `QuantumScript`, as it is no longer necessary. There is
   still a `_qfunc_output` property on `QNode` instances.
   [(#4651)](https://github.com/PennyLaneAI/pennylane/pull/4651)
@@ -198,6 +203,7 @@
   of the computed qubit operator, if imaginary components are smaller than a threshold. 
   [(#4639)](https://github.com/PennyLaneAI/pennylane/pull/4639)
 
+
 <h3>Breaking changes 💔</h3>
 
 * The device test suite now converts device kwargs to integers or floats if they can be converted to integers or floats.
diff --git a/pennylane/transforms/decompositions/single_qubit_unitary.py b/pennylane/transforms/decompositions/single_qubit_unitary.py
index b19bbc8daa9..62969e36b7a 100644
--- a/pennylane/transforms/decompositions/single_qubit_unitary.py
+++ b/pennylane/transforms/decompositions/single_qubit_unitary.py
@@ -15,7 +15,7 @@
 operations into elementary gates.
 """
 
-import numpy
+import numpy as np
 
 import pennylane as qml
 from pennylane import math
@@ -40,9 +40,9 @@ def _convert_to_su2(U, return_global_phase=False):
     # Compute the determinants
     U = qml.math.cast(U, "complex128")
     dets = math.linalg.det(U)
+    exp_angles = math.angle(dets) / 2
+    U_SU2 = math.cast_like(U, dets) * math.exp(-1j * math.cast_like(exp_angles, 1j))[:, None, None]
 
-    exp_angles = math.cast_like(math.angle(dets), 1j) / 2
-    U_SU2 = math.cast_like(U, dets) * math.exp(-1j * exp_angles)[:, None, None]
     return (U_SU2, exp_angles) if return_global_phase else U_SU2
 
 
@@ -186,9 +186,9 @@ def _zyz_decomposition(U, wire, return_global_phase=False):
 
     phis, thetas, omegas, alphas = map(math.squeeze, [phis, thetas, omegas, alphas])
 
-    phis = phis % (4 * numpy.pi)
-    thetas = thetas % (4 * numpy.pi)
-    omegas = omegas % (4 * numpy.pi)
+    phis = phis % (4 * np.pi)
+    thetas = thetas % (4 * np.pi)
+    omegas = omegas % (4 * np.pi)
 
     operations = [qml.RZ(phis, wire), qml.RY(thetas, wire), qml.RZ(omegas, wire)]
     if return_global_phase:
@@ -250,9 +250,9 @@ def _xyx_decomposition(U, wire, return_global_phase=False):
 
     phis, thetas, lams, gammas = map(math.squeeze, [phis, thetas, lams, gammas])
 
-    phis = phis % (4 * numpy.pi)
-    thetas = thetas % (4 * numpy.pi)
-    lams = lams % (4 * numpy.pi)
+    phis = phis % (4 * np.pi)
+    thetas = thetas % (4 * np.pi)
+    lams = lams % (4 * np.pi)
 
     operations = [qml.RX(lams, wire), qml.RY(thetas, wire), qml.RX(phis, wire)]
     if return_global_phase:
@@ -316,9 +316,9 @@ def _zxz_decomposition(U, wire, return_global_phase=False):
 
     phis, thetas, psis, alphas = map(math.squeeze, [phis, thetas, psis, alphas])
 
-    phis = phis % (4 * numpy.pi)
-    thetas = thetas % (4 * numpy.pi)
-    psis = psis % (4 * numpy.pi)
+    phis = phis % (4 * np.pi)
+    thetas = thetas % (4 * np.pi)
+    psis = psis % (4 * np.pi)
 
     # Return gates in the order they will be applied on the qubit
     operations = [qml.RZ(psis, wire), qml.RX(thetas, wire), qml.RZ(phis, wire)]
diff --git a/tests/transforms/test_decompositions.py b/tests/transforms/test_decompositions.py
index ccaaafd315f..b4bdda9688d 100644
--- a/tests/transforms/test_decompositions.py
+++ b/tests/transforms/test_decompositions.py
@@ -38,16 +38,16 @@
 typeof_gates_zyz = (qml.RZ, qml.RY, qml.RZ, qml.GlobalPhase)
 single_qubit_decomps_zyz = [
     (I, typeof_gates_zyz, [0.0, 0.0, 0.0, 0]),
-    (Z, typeof_gates_zyz, [np.pi / 2, 0.0, np.pi / 2, (-1.5707963267948966 + 0j)]),
+    (Z, typeof_gates_zyz, [np.pi / 2, 0.0, np.pi / 2, -1.5707963267948966]),
     (
         S,
         typeof_gates_zyz,
-        [np.pi / 4, 0.0, np.pi / 4, (-0.7853981633974483 - 1.6780315470477092e-09j)],
+        [np.pi / 4, 0.0, np.pi / 4, -0.7853981633974483],
     ),
     (
         T,
         typeof_gates_zyz,
-        [np.pi / 8, 0.0, np.pi / 8, (-0.39269908047469393 - 3.225207101387184e-09j)],
+        [np.pi / 8, 0.0, np.pi / 8, -0.39269908047469393],
     ),
     (qml.RZ(0.3, wires=0).matrix(), typeof_gates_zyz, [0.15, 0.0, 0.15, 0]),
     (
@@ -70,8 +70,8 @@
         typeof_gates_zyz,
         [12.382273469673908, np.pi, 0.18409714468526372, 0],
     ),
-    (H, typeof_gates_zyz, [np.pi, np.pi / 2, 0.0, (-1.5707963267948966 + 0j)]),
-    (X, typeof_gates_zyz, [np.pi / 2, np.pi, 10.995574287564276, (-1.5707963267948966 + 0j)]),
+    (H, typeof_gates_zyz, [np.pi, np.pi / 2, 0.0, -1.5707963267948966]),
+    (X, typeof_gates_zyz, [np.pi / 2, np.pi, 10.995574287564276, -1.5707963267948966]),
     (
         np.exp(1j * 0.02) * qml.Rot(-1.0, 2.0, -3.0, wires=0).matrix(),
         typeof_gates_zyz,
@@ -79,7 +79,7 @@
             11.566370614359172,
             2.0,
             9.566370614359172,
-            (-0.020000000000000042 - 2.122325752640375e-17j),
+            -0.020000000000000042,
         ],
     ),
     # Add two instances of broadcasted unitaries, one coming from RZ and another from Rot
@@ -196,14 +196,14 @@ def test_zyz_decomposition_jax(self, U, expected_gates, expected_params):
             10.845351366405708,
             1.3974974118006183,
             0.45246583660683803,
-            (1.1759220332464762 - 4.163336342344337e-17j),
+            1.1759220332464762,
         ),
     ),
     # Try a few specific special unitaries
     (I, typeof_gates_xyx, [0, 0, 0, 0]),  # This triggers the if conditional trivially
-    (X, typeof_gates_xyx, [4.71238898038469, 0.0, 10.995574287564276, (-1.5707963267948966 + 0j)]),
-    (Y, typeof_gates_xyx, [1 / 2 * np.pi, np.pi, 1 / 2 * np.pi, (-1.5707963267948966 + 0j)]),
-    (Z, typeof_gates_xyx, [10.995574287564276, np.pi, 1 / 2 * np.pi, (-1.5707963267948966 + 0j)]),
+    (X, typeof_gates_xyx, [4.71238898038469, 0.0, 10.995574287564276, -1.5707963267948966]),
+    (Y, typeof_gates_xyx, [1 / 2 * np.pi, np.pi, 1 / 2 * np.pi, -1.5707963267948966]),
+    (Z, typeof_gates_xyx, [10.995574287564276, np.pi, 1 / 2 * np.pi, -1.5707963267948966]),
     # Add two instances of broadcasted unitaries, one coming from RZ and another from Rot
     (
         qml.QubitUnitary(qml.RZ.compute_matrix(np.array([np.pi, np.pi / 2])), wires=0).matrix(),
@@ -311,16 +311,16 @@ def test_xyx_decomposition_jax(self, U, expected_gates, expected_params):
 typeof_gates_zxz = (qml.RZ, qml.RX, qml.RZ, qml.GlobalPhase)
 single_qubit_decomps_zxz = [
     (I, typeof_gates_zxz, [0.0, 0.0, 0.0, 0]),
-    (Z, typeof_gates_zxz, [np.pi / 2, 0.0, np.pi / 2, (-1.5707963267948966 + 0j)]),
+    (Z, typeof_gates_zxz, [np.pi / 2, 0.0, np.pi / 2, -1.5707963267948966]),
     (
         S,
         typeof_gates_zxz,
-        [np.pi / 4, 0.0, np.pi / 4, (-0.7853981633974483 - 1.6780315470477092e-09j)],
+        [np.pi / 4, 0.0, np.pi / 4, -0.7853981633974483],
     ),
     (
         T,
         typeof_gates_zxz,
-        [np.pi / 8, 0.0, np.pi / 8, (-0.39269908047469393 - 3.225207101387184e-09j)],
+        [np.pi / 8, 0.0, np.pi / 8, -0.39269908047469393],
     ),
     (qml.RZ(0.3, wires=0).matrix(), typeof_gates_zxz, [0.15, 0.0, 0.15, 0]),
     (
@@ -484,7 +484,7 @@ def test_zxz_decomposition_jax(self, U, expected_gates, expected_params):
             10.845351366405708,
             1.3974974118006183,
             0.45246583660683803,
-            1.1759220332464762 - 4.163336342344337e-17j,
+            1.1759220332464762,
         ),
     ),
     (