diff --git a/qiskit/synthesis/discrete_basis/generate_basis_approximations.py b/qiskit/synthesis/discrete_basis/generate_basis_approximations.py
index 672d0eb9e8ef..da9708c24559 100644
--- a/qiskit/synthesis/discrete_basis/generate_basis_approximations.py
+++ b/qiskit/synthesis/discrete_basis/generate_basis_approximations.py
@@ -156,7 +156,7 @@ def generate_basic_approximations(
         data = {}
         for sequence in sequences:
             gatestring = sequence.name
-            data[gatestring] = sequence.product
+            data[gatestring] = (sequence.product, sequence.global_phase)
 
         np.save(filename, data)
 
diff --git a/qiskit/synthesis/discrete_basis/solovay_kitaev.py b/qiskit/synthesis/discrete_basis/solovay_kitaev.py
index 2c8df5bd1b6d..f367f6c0f0b5 100644
--- a/qiskit/synthesis/discrete_basis/solovay_kitaev.py
+++ b/qiskit/synthesis/discrete_basis/solovay_kitaev.py
@@ -51,14 +51,19 @@ def __init__(
 
         self.basic_approximations = self.load_basic_approximations(basic_approximations)
 
-    def load_basic_approximations(self, data: list | str | dict) -> list[GateSequence]:
+    @staticmethod
+    def load_basic_approximations(data: list | str | dict) -> list[GateSequence]:
         """Load basic approximations.
 
         Args:
             data: If a string, specifies the path to the file from where to load the data.
-                If a dictionary, directly specifies the decompositions as ``{gates: matrix}``.
-                There ``gates`` are the names of the gates producing the SO(3) matrix ``matrix``,
-                e.g. ``{"h t": np.array([[0, 0.7071, -0.7071], [0, -0.7071, -0.7071], [-1, 0, 0]]}``.
+                If a dictionary, directly specifies the decompositions as ``{gates: matrix}``
+                or ``{gates: (matrix, global_phase)}``. There, ``gates`` are the names of the gates
+                producing the SO(3) matrix ``matrix``, e.g.
+                ``{"h t": np.array([[0, 0.7071, -0.7071], [0, -0.7071, -0.7071], [-1, 0, 0]]}``
+                and the ``global_phase`` can be given to account for a global phase difference
+                between the U(2) matrix of the quantum gates and the stored SO(3) matrix.
+                If not given, the ``global_phase`` will be assumed to be 0.
 
         Returns:
             A list of basic approximations as type ``GateSequence``.
@@ -72,13 +77,20 @@ def load_basic_approximations(self, data: list | str | dict) -> list[GateSequenc
 
         # if a file, load the dictionary
         if isinstance(data, str):
-            data = np.load(data, allow_pickle=True)
+            data = np.load(data, allow_pickle=True).item()
 
         sequences = []
-        for gatestring, matrix in data.items():
+        for gatestring, matrix_and_phase in data.items():
+            if isinstance(matrix_and_phase, tuple):
+                matrix, global_phase = matrix_and_phase
+            else:
+                matrix, global_phase = matrix_and_phase, 0
+
             sequence = GateSequence()
             sequence.gates = [_1q_gates[element] for element in gatestring.split()]
+            sequence.labels = [gate.name for gate in sequence.gates]
             sequence.product = np.asarray(matrix)
+            sequence.global_phase = global_phase
             sequences.append(sequence)
 
         return sequences
diff --git a/releasenotes/notes/fix-sk-load-from-file-02c6eabbbd7fcda3.yaml b/releasenotes/notes/fix-sk-load-from-file-02c6eabbbd7fcda3.yaml
new file mode 100644
index 000000000000..d995af06bccb
--- /dev/null
+++ b/releasenotes/notes/fix-sk-load-from-file-02c6eabbbd7fcda3.yaml
@@ -0,0 +1,10 @@
+---
+fixes:
+  - |
+    Fix the :class:`.SolovayKitaev` transpiler pass when loading basic
+    approximations from an exising ``.npy`` file. Previously, loading
+    a stored approximation which allowed for further reductions (e.g. due
+    to gate cancellations) could cause a runtime failure.
+    Additionally, the global phase difference of the U(2) gate product
+    and SO(3) representation was lost during a save-reload procedure.
+    Fixes `Qiskit/qiskit#12576 <https://github.com/Qiskit/qiskit/issues/12576>`_.
diff --git a/test/python/transpiler/test_solovay_kitaev.py b/test/python/transpiler/test_solovay_kitaev.py
index e15a080f6f03..62b811c8e3bd 100644
--- a/test/python/transpiler/test_solovay_kitaev.py
+++ b/test/python/transpiler/test_solovay_kitaev.py
@@ -12,8 +12,10 @@
 
 """Test the Solovay Kitaev transpilation pass."""
 
+import os
 import unittest
 import math
+import tempfile
 import numpy as np
 import scipy
 
@@ -230,6 +232,35 @@ def test_u_gates_work(self):
         included_gates = set(discretized.count_ops().keys())
         self.assertEqual(set(basis_gates), included_gates)
 
+    def test_load_from_file(self):
+        """Test loading basic approximations from a file works.
+
+        Regression test of Qiskit/qiskit#12576.
+        """
+        filename = "approximations.npy"
+
+        with tempfile.TemporaryDirectory() as tmp_dir:
+            fullpath = os.path.join(tmp_dir, filename)
+
+            # dump approximations to file
+            generate_basic_approximations(basis_gates=["h", "s", "sdg"], depth=3, filename=fullpath)
+
+            # circuit to decompose and reference decomp
+            circuit = QuantumCircuit(1)
+            circuit.rx(0.8, 0)
+
+            reference = QuantumCircuit(1, global_phase=3 * np.pi / 4)
+            reference.h(0)
+            reference.s(0)
+            reference.h(0)
+
+            # load the decomp and compare to reference
+            skd = SolovayKitaev(basic_approximations=fullpath)
+            # skd = SolovayKitaev(basic_approximations=filename)
+            discretized = skd(circuit)
+
+        self.assertEqual(discretized, reference)
+
 
 @ddt
 class TestGateSequence(QiskitTestCase):