Repository: Qiskit/qiskit (python)

Issue:
Solovay Kitaev Decomposition doesn't correctly load saved approximations from file
### Environment

- **Qiskit version**: 1.1.0
- **Python version**: 3.10.11
- **Operating system**: Windows 11 23H2
- **Numpy version**: 1.26.4

### What is happening?

After specifying a file to save the approximations of `generate_basic_approximations` into, attempting to load this into a new `SolovayKitaevDecomposition` gives a runtime error.

This is somewhat easily fixable by replacing [this line](https://github.com/Qiskit/qiskit/blob/03aeee81fefb2f915a82ceee37d64971eadbf6f0/qiskit/synthesis/discrete_basis/solovay_kitaev.py#L80) with the following temporary fix:
```py
        for gatestring, matrix in data.flat[0].items():
```
After this, the `SolovayKitaevDecomposition` object can be instantiated, but still contains a separate logic error. Some inputs can cause the circuits generated by this object to have different global phases than the ones generated by a `SolovayKitaevDecomposition` which did not load these approximations from a file (see example).

### How can we reproduce the issue?

## Minimal Example
```py
from qiskit.synthesis.discrete_basis.solovay_kitaev import (
    SolovayKitaevDecomposition,
)
from qiskit.synthesis.discrete_basis import generate_basic_approximations

import numpy as np
from qiskit.synthesis.discrete_basis.gate_sequence import (
    _convert_so3_to_su2,
)

generate_basic_approximations(["t", "tdg", "h"], 10, "test.npy")
sk = SolovayKitaevDecomposition()
sk_file = SolovayKitaevDecomposition("test.npy")
# first issue occurs on line above

rotation = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])

qc = sk.run(_convert_so3_to_su2(rotation), 2)
qc_file = sk_file.run(_convert_so3_to_su2(rotation), 2)

print(f"{qc.global_phase = }")
print(f"{qc_file.global_phase = }")

print(f"qc =\n{qc.draw()}")
print(f"qc_file =\n{qc_file.draw()}")
```
## Output
### Runtime Error on loading file
```
Traceback (most recent call last):
  File "[omitted]\saved_approximations.py", line 13, in <module>
    sk_file = SolovayKitaevDecomposition("test.npy")    
  File "[omitted]\env\lib\site-packages\qiskit\synthesis\discrete_basis\solovay_kitaev.py", line 54, in __init__
    self.basic_approximations = self.load_basic_approximations(basic_approximations)
  File "[omitted]\env\lib\site-packages\qiskit\synthesis\discrete_basis\solovay_kitaev.py", line 80, in load_basic_approximations
    for gatestring, matrix in data.items():
AttributeError: 'numpy.ndarray' object has no attribute 'items'. Did you mean: 'item'?
```
### After Runtime Error Fix
```
qc.global_phase = 3.141592653589793
qc_file.global_phase = 0.0
qc =
global phase: π
   ┌───┐┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ T ├┤ T ├┤ T ├┤ T ├
   └───┘└───┘└───┘└───┘└───┘
qc_file =
   ┌───┐┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ T ├┤ T ├┤ T ├┤ T ├
   └───┘└───┘└───┘└───┘└───┘ 
```

### What should happen?

The circuit created using the `SolovayKitaevDecomposition` created with saved approximations (`qc_file`) should behave the same as the one without (`qc`) - they should both give the same result gates and global phase, as well as run without runtime errors.

### Any suggestions?

- The temporary fix for the `AttributeError` issue should probably not be used (or at least checked before used). On my system, at least, the value of `data` at line 79 in `solovay_kitaev.py` is a 0 dimensional `np.ndarray` containing the object loaded from the file. I am not sure whether this is consistent on other systems or numpy versions.
  - Also, `data.size` is not checked in this one-line solution, so there could be some sort of security vulnerability if this is used.
- It seems that only the gates and product are saved into the file ([source code here](https://github.com/Qiskit/qiskit/blob/03aeee81fefb2f915a82ceee37d64971eadbf6f0/qiskit/synthesis/discrete_basis/generate_basis_approximations.py#L157)). It is also mentioned explicitly in the [constructor docstring of `SolovayKitaevDecomposition`](https://github.com/Qiskit/qiskit/blob/03aeee81fefb2f915a82ceee37d64971eadbf6f0/qiskit/synthesis/discrete_basis/solovay_kitaev.py#L46), so I'm unsure if this is intentional. If not and this needs to be changed, this documentation should also be updated.

Pull request:
Fix some bugs in loading Solovay Kitaev decompositions (backport #12579)


### Summary

Fixes #12576 and some other small bugs in storing & loading basic approximations for the Solovay Kitaev algorithm.

### Details and comments

* as #12576 mentioned, there is a missing `.item()` unwrapping when loading the dictionary from the `npy` file
* reloading the gate sequence also missed the correct `labels` and `global_phase`
* a test is now checking all of the above works correctly
<hr>This is an automatic backport of pull request #12579 done by [Mergify](https://mergify.com).

Patch:
diff --git a/qiskit/synthesis/discrete_basis/generate_basis_approximations.py b/qiskit/synthesis/discrete_basis/generate_basis_approximations.py
index 07139b223..73ee70446 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 c0cfbd81c..f6af0d667 100644
--- a/qiskit/synthesis/discrete_basis/solovay_kitaev.py
+++ b/qiskit/synthesis/discrete_basis/solovay_kitaev.py
@@ -51,14 +51,19 @@ class SolovayKitaevDecomposition:
 
         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 @@ class SolovayKitaevDecomposition:
 
         # 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 000000000..d995af06b
--- /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>`_.


Test patch:
diff --git a/test/python/transpiler/test_solovay_kitaev.py b/test/python/transpiler/test_solovay_kitaev.py
index e15a080f6..62b811c8e 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 @@ class TestSolovayKitaev(QiskitTestCase):
         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):