diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 0f65c68..55fd90b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -16,10 +16,11 @@ repos: args: # compatibility with black # E203 whitespace before ':' + # E704 multiple statements on one line (def) # E741 ambiguous variable name 'I' # W503 line break before binary operator - "--max-line-length=88" - - "--ignore=E203, E741, W503" + - "--ignore=E203, E704, E741, W503" - repo: https://github.com/pre-commit/mirrors-mypy rev: 'v1.15.0' hooks: @@ -31,7 +32,8 @@ repos: name: pytest entry: bash -c 'PYTHONPATH=./ .venv/bin/pytest tests' language: system - types: [python] + pass_filenames: false + always_run: true - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.4.0 hooks: diff --git a/.vscode/settings.json b/.vscode/settings.json index cbe7cee..8506f48 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -4,65 +4,92 @@ }, "editor.formatOnSave": true, "cSpell.words": [ + "addopts", + "allclose", "Amano", + "capsys", + "choi", "CINV", - "Kazuyuki", - "Kliuchnikov", - "Maslov", - "Mosca", - "Matsumoto", - "Nobuyuki", - "ODGP", - "PYTHONPATH", - "Selinger", - "Shuntaro", - "TCONJ", - "TDGP", - "Vadym", - "Yamamoto", - "Yoshioka", - "addopts", "cnot", + "cnots", + "cvxpy", "decomp", "denomexp", "dloop", "domega", "droottwo", "dtimeout", + "ECOS", + "eigh", "eigsy", + "eigval", "facs", "figsize", "floop", "floorlog", "floorsqrt", "ftimeout", + "FWHT", + "gptm", "gridsynth", + "GUROBI", "isort", + "Kazuyuki", + "Kliuchnikov", + "kron", + "kwargs", "limegreen", + "linalg", + "Maslov", "matplotlib", + "Matsumoto", + "Mosca", "mpmath", "mymath", "mypy", + "ndarray", "newsynth", + "njit", + "Nobuyuki", + "nonneg", + "numba", "numpy", + "ODGP", "orangered", "prec", + "prefactors", + "probs", "pygridsynth", "pypi", "pytest", + "PYTHONPATH", "qubit", "qubits", + "qulacs", + "randn", "rounddiv", + "scipy", "selfassociate", "selfcoprime", + "Selinger", + "semidefinite", "setuptools", "showgraph", + "Shuntaro", + "TCONJ", + "TDGP", "testpaths", + "triu", "unitaries", + "Vadym", + "vecs", "venv", + "Weyl", "workdps", "xlabel", + "Yamamoto", "ylabel", + "Yoshioka", "zomega", "zroottwo", ], diff --git a/README.md b/README.md index 8774f7e..ee4b39c 100644 --- a/README.md +++ b/README.md @@ -69,7 +69,7 @@ pygridsynth [options] - `--dloop`, `-dl`: Maximum number of failed integer factoring attempts allowed during Diophantine equation solving (default: `10`). - `--floop`, `-fl`: Maximum number of failed integer factoring attempts allowed during the factoring process (default: `10`). - `--seed`: Random seed for deterministic results.(default: `0`) -- `--verbose`, `-v`: Enables detailed output. +- `--verbose`, `-v`: Verbosity level (0=silent, 1=basic, 2=detailed, 3=debug).(default: `0`) - `--time`, `-t`: Measures the execution time. - `--showgraph`, `-g`: Displays the decomposition result as a graph. - `--up-to-phase`, `-ph`: Approximates up to a phase. @@ -138,6 +138,157 @@ gates = gridsynth_gates(theta=theta, epsilon=epsilon) print(gates) ``` +### Multi-Qubit Unitary Approximation + +`pygridsynth` provides functionality for approximating multi-qubit unitary matrices using the Clifford+T gate set. This is useful for synthesizing quantum circuits that implement arbitrary multi-qubit unitaries. + +**Basic usage:** + +```python +import mpmath + +from pygridsynth.multi_qubit_unitary_approximation import ( + approximate_multi_qubit_unitary, +) + +# Define a target unitary matrix (example: 2-qubit identity) +num_qubits = 2 +U = mpmath.eye(2**num_qubits) # 4x4 identity matrix +epsilon = "1e-10" + +# Approximate the unitary +circuit, U_approx = approximate_multi_qubit_unitary(U, num_qubits, epsilon) + +print(f"Circuit length: {len(circuit)}") +print(f"Circuit: {str(circuit)}") +``` + +**Using with random unitary:** + +```python +from pygridsynth.multi_qubit_unitary_approximation import ( + approximate_multi_qubit_unitary, +) +from pygridsynth.mymath import random_su + +# Generate a random SU(2^n) unitary +num_qubits = 2 +U = random_su(num_qubits) +epsilon = "1e-10" + +# Approximate with high precision +circuit, U_approx = approximate_multi_qubit_unitary(U, num_qubits, epsilon) +``` + +**Returning DOmegaMatrix:** + +```python +from pygridsynth.multi_qubit_unitary_approximation import ( + approximate_multi_qubit_unitary, +) +from pygridsynth.mymath import random_su + +# Generate a random SU(2^n) unitary +num_qubits = 2 +U = random_su(num_qubits) +epsilon = "1e-10" + +# Return DOmegaMatrix instead of mpmath.matrix for more efficient representation +circuit, U_domega = approximate_multi_qubit_unitary( + U, num_qubits, epsilon, return_domega_matrix=True +) + +# Convert to complex matrix if needed +U_complex = U_domega.to_complex_matrix +``` + +**Parameters:** + +- `U`: Target unitary matrix (`mpmath.matrix`) +- `num_qubits`: Number of qubits +- `epsilon`: Error tolerance (can be `str`, `float`, or `mpmath.mpf`) +- `return_domega_matrix`: If `True`, returns `DOmegaMatrix`; if `False`, returns `mpmath.matrix` (default: `False`) +- `scale_epsilon`: Whether to scale epsilon based on the number of qubits (default: `True`) +- `cfg`: Optional `GridsynthConfig` object for advanced configuration +- `**kwargs`: Additional configuration options (ignored if `cfg` is provided) + +**Returns:** + +A tuple of `(circuit, U_approx)`: +- `circuit`: `QuantumCircuit` object representing the Clifford+T decomposition +- `U_approx`: Approximated unitary matrix (`mpmath.matrix` or `DOmegaMatrix` depending on `return_domega_matrix`) + +### Mixed Unitary Synthesis + +`pygridsynth` also provides functionality for mixed unitary synthesis, which approximates a target unitary by mixing multiple perturbed unitaries. This is useful for reducing the number of T-gates in quantum circuits. + +The library provides two versions: `mixed_synthesis_parallel` (for parallel execution) and `mixed_synthesis_sequential` (for sequential execution). + +**Basic usage with mpmath.matrix:** + +```python +from pygridsynth.mixed_synthesis import mixed_synthesis_sequential +from pygridsynth.mymath import diamond_norm_error_from_choi, random_su + +# Generate a random SU(2^n) unitary matrix +num_qubits = 2 +unitary = random_su(num_qubits) + +# Parameters +eps = 1e-4 # Error tolerance +M = 64 # Number of Hermitian operators for perturbation +seed = 123 # Random seed for reproducibility + +# Compute mixed synthesis (sequential version) +result = mixed_synthesis_sequential(unitary, num_qubits, eps, M, seed=seed) + +if result is not None: + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + print(f"Number of circuits: {len(circuit_list)}") + print(f"Mixing probabilities: {probs_gptm}") + error = diamond_norm_error_from_choi(u_choi, u_choi_opt, eps, mixed_synthesis=True) + print(f"error: {error}") +``` + +**Using parallel version:** + +```python +import mpmath + +from pygridsynth.mixed_synthesis import mixed_synthesis_parallel + +# Generate a random SU(2^n) unitary matrix +num_qubits = 2 +unitary = mpmath.eye(2**num_qubits) + +# Parameters +eps = 1e-4 # Error tolerance +M = 64 # Number of Hermitian operators for perturbation +seed = 123 # Random seed for reproducibility + +# For faster computation with multiple cores +result = mixed_synthesis_parallel(unitary, num_qubits, eps, M, seed=seed) +``` + +**Parameters:** + +- `unitary`: Target unitary matrix (`mpmath.matrix` or `numpy.ndarray`) +- `num_qubits`: Number of qubits +- `eps`: Error tolerance parameter +- `M`: Number of Hermitian operators for perturbation +- `seed`: Random seed for reproducibility (default: `123`) +- `dps`: Decimal precision (default: `-1` for auto-calculation) + +**Returns:** + +A tuple of `(circuit_list, eu_np_list, probs_gptm, u_choi_opt)` or `None` on failure: +- `circuit_list`: List of `QuantumCircuit` objects for perturbed unitaries +- `eu_np_list`: List of approximated unitary matrices (numpy arrays) +- `probs_gptm`: Array of mixing probabilities +- `u_choi_opt`: Optimal mixed Choi matrix + +**Note:** The parallel version (`mixed_synthesis_parallel`) uses multiprocessing and may be faster for large `M` values, while the sequential version (`mixed_synthesis_sequential`) is more suitable for debugging or when parallel execution is not desired. + ## Contributing @@ -149,9 +300,15 @@ This project is licensed under the MIT License. ## References -- Brett Giles and Peter Selinger. Remarks on Matsumoto and Amano's normal form for single-qubit Clifford+T operators, 2019. -- Ken Matsumoto and Kazuyuki Amano. Representation of Quantum Circuits with Clifford and π/8 Gates, 2008. -- Neil J. Ross and Peter Selinger. Optimal ancilla-free Clifford+T approximation of z-rotations, 2016. -- Peter Selinger. Efficient Clifford+T approximation of single-qubit operators, 2014. -- Peter Selinger and Neil J. Ross. Exact and approximate synthesis of quantum circuits. https://www.mathstat.dal.ca/~selinger/newsynth/, 2018. -- Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca. Fast and efficient exact synthesis of single qubit unitaries generated by Clifford and T gates, 2013. +- Vadym Kliuchnikov, Kristin Lauter, Romy Minko, Adam Paetznick, Christophe Petit. "Shorter quantum circuits via single-qubit gate approximation." Quantum 7 (2023): 1208. DOI: 10.22331/q-2023-12-18-1208. +- Peter Selinger. "Efficient Clifford+T approximation of single-qubit operators." Quantum Info. Comput. 15, no. 1-2 (2015): 159-180. +- Neil J. Ross and Peter Selinger. "Optimal ancilla-free Clifford+T approximation of z-rotations." Quantum Info. Comput. 16, no. 11-12 (2016): 901-953. +- Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca. "Fast and efficient exact synthesis of single-qubit unitaries generated by Clifford and T gates." Quantum Info. Comput. 13, no. 7-8 (2013): 607-630. +- Ken Matsumoto and Kazuyuki Amano. "Representation of Quantum Circuits with Clifford and π/8 Gates." arXiv: Quantum Physics (2008). URL: https://api.semanticscholar.org/CorpusID:17327793. +- Brett Gordon Giles and Peter Selinger. "Remarks on Matsumoto and Amano's normal form for single-qubit Clifford+T operators." ArXiv abs/1312.6584 (2013). URL: https://api.semanticscholar.org/CorpusID:10077777. +- Vivek V. Shende, Igor L. Markov, and Stephen S. Bullock. "Minimal universal two-qubit controlled-NOT-based circuits." Phys. Rev. A 69, no. 6 (2004): 062321. DOI: 10.1103/PhysRevA.69.062321. +- Vivek V. Shende, Igor L. Markov, and Stephen S. Bullock. "Finding Small Two-Qubit Circuits." Proceedings of SPIE - The International Society for Optical Engineering (2004). DOI: 10.1117/12.542381. +- Vivek V. Shende, Igor L. Markov, and Stephen S. Bullock. "Smaller Two-Qubit Circuits for Quantum Communication and Computation." In Proceedings of the Conference on Design, Automation and Test in Europe - Volume 2, p. 20980. IEEE Computer Society, 2004. +- Korbinian Kottmann. "Two-qubit Synthesis." (2025). URL: https://pennylane.ai/compilation/two-qubit-synthesis. +- Anna M. Krol and Zaid Al-Ars. "Beyond quantum Shannon decomposition: Circuit construction for n-qubit gates based on block-ZXZ decomposition." Physical Review Applied 22, no. 3 (2024): 034019. DOI: 10.1103/PhysRevApplied.22.034019. +- Peter Selinger and Neil J. Ross. "Exact and approximate synthesis of quantum circuits." (2018). URL: https://www.mathstat.dal.ca/~selinger/newsynth/. diff --git a/dist/pygridsynth-1.2.0-py3-none-any.whl b/dist/pygridsynth-1.2.0-py3-none-any.whl deleted file mode 100644 index 8c5b326..0000000 Binary files a/dist/pygridsynth-1.2.0-py3-none-any.whl and /dev/null differ diff --git a/dist/pygridsynth-1.2.0.tar.gz b/dist/pygridsynth-1.2.0.tar.gz deleted file mode 100644 index 6c1bd89..0000000 Binary files a/dist/pygridsynth-1.2.0.tar.gz and /dev/null differ diff --git a/examples/mixed_synthesis_parallel.py b/examples/mixed_synthesis_parallel.py new file mode 100644 index 0000000..ed1c2d2 --- /dev/null +++ b/examples/mixed_synthesis_parallel.py @@ -0,0 +1,15 @@ +import mpmath + +from pygridsynth.mixed_synthesis import mixed_synthesis_parallel + +# Generate a random SU(2^n) unitary matrix +num_qubits = 2 +unitary = mpmath.eye(2**num_qubits) + +# Parameters +eps = 1e-4 # Error tolerance +M = 64 # Number of Hermitian operators for perturbation +seed = 123 # Random seed for reproducibility + +# For faster computation with multiple cores +result = mixed_synthesis_parallel(unitary, num_qubits, eps, M, seed=seed) diff --git a/examples/mixed_synthesis_sequential.py b/examples/mixed_synthesis_sequential.py new file mode 100644 index 0000000..0b9b1a2 --- /dev/null +++ b/examples/mixed_synthesis_sequential.py @@ -0,0 +1,21 @@ +from pygridsynth.mixed_synthesis import mixed_synthesis_sequential +from pygridsynth.mymath import diamond_norm_error_from_choi, random_su + +# Generate a random SU(2^n) unitary matrix +num_qubits = 2 +unitary = random_su(num_qubits) + +# Parameters +eps = 1e-4 # Error tolerance +M = 64 # Number of Hermitian operators for perturbation +seed = 123 # Random seed for reproducibility + +# Compute mixed synthesis (sequential version) +result = mixed_synthesis_sequential(unitary, num_qubits, eps, M, seed=seed) + +if result is not None: + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + print(f"Number of circuits: {len(circuit_list)}") + print(f"Mixing probabilities: {probs_gptm}") + error = diamond_norm_error_from_choi(u_choi, u_choi_opt, eps, mixed_synthesis=True) + print(f"error: {error}") diff --git a/examples/multi_qubit_basic.py b/examples/multi_qubit_basic.py new file mode 100644 index 0000000..765ec48 --- /dev/null +++ b/examples/multi_qubit_basic.py @@ -0,0 +1,16 @@ +import mpmath + +from pygridsynth.multi_qubit_unitary_approximation import ( + approximate_multi_qubit_unitary, +) + +# Define a target unitary matrix (example: 2-qubit identity) +num_qubits = 2 +U = mpmath.eye(2**num_qubits) # 4x4 identity matrix +epsilon = "1e-10" + +# Approximate the unitary +circuit, U_approx = approximate_multi_qubit_unitary(U, num_qubits, epsilon) + +print(f"Circuit length: {len(circuit)}") +print(f"Circuit: {str(circuit)}") diff --git a/examples/multi_qubit_domega.py b/examples/multi_qubit_domega.py new file mode 100644 index 0000000..288c8db --- /dev/null +++ b/examples/multi_qubit_domega.py @@ -0,0 +1,17 @@ +from pygridsynth.multi_qubit_unitary_approximation import ( + approximate_multi_qubit_unitary, +) +from pygridsynth.mymath import random_su + +# Generate a random SU(2^n) unitary +num_qubits = 2 +U = random_su(num_qubits) +epsilon = "1e-10" + +# Return DOmegaMatrix instead of mpmath.matrix for more efficient representation +circuit, U_domega = approximate_multi_qubit_unitary( + U, num_qubits, epsilon, return_domega_matrix=True +) + +# Convert to complex matrix if needed +U_complex = U_domega.to_complex_matrix diff --git a/examples/multi_qubit_random.py b/examples/multi_qubit_random.py new file mode 100644 index 0000000..e1f18d1 --- /dev/null +++ b/examples/multi_qubit_random.py @@ -0,0 +1,12 @@ +from pygridsynth.multi_qubit_unitary_approximation import ( + approximate_multi_qubit_unitary, +) +from pygridsynth.mymath import random_su + +# Generate a random SU(2^n) unitary +num_qubits = 2 +U = random_su(num_qubits) +epsilon = "1e-10" + +# Approximate with high precision +circuit, U_approx = approximate_multi_qubit_unitary(U, num_qubits, epsilon) diff --git a/pygridsynth/__main__.py b/pygridsynth/__main__.py index 8ac78a9..9ae637f 100644 --- a/pygridsynth/__main__.py +++ b/pygridsynth/__main__.py @@ -1,4 +1,4 @@ from .cli import main if __name__ == "__main__": - print(main()) + main() diff --git a/pygridsynth/cli.py b/pygridsynth/cli.py index e6fd965..77c216d 100644 --- a/pygridsynth/cli.py +++ b/pygridsynth/cli.py @@ -15,10 +15,11 @@ "allowed during the factoring process" ), "seed": "Random seed for deterministic results", + "verbose": "Verbosity level (0=silent, 1=basic, 2=detailed, 3=debug)", } -def main() -> str: +def main() -> None: parser = argparse.ArgumentParser() parser.add_argument("theta", type=str) @@ -29,7 +30,7 @@ def main() -> str: parser.add_argument("--dloop", "-dl", type=int, help=helps["dl"]) parser.add_argument("--floop", "-fl", type=int, help=helps["fl"]) parser.add_argument("--seed", type=int, help=helps["seed"]) - parser.add_argument("--verbose", "-v", action="store_true") + parser.add_argument("--verbose", "-v", type=int) parser.add_argument("--time", "-t", action="store_true") parser.add_argument("--showgraph", "-g", action="store_true") parser.add_argument("--up-to-phase", "-ph", action="store_true") @@ -45,7 +46,7 @@ def main() -> str: if args.dloop is not None: cfg_args["dloop"] = args.dloop if args.floop is not None: - cfg_args["floop"] = args + cfg_args["floop"] = args.floop if args.seed is not None: cfg_args["seed"] = args.seed if args.verbose: @@ -60,4 +61,5 @@ def main() -> str: cfg = GridsynthConfig(**cfg_args) gates = gridsynth_gates(theta=args.theta, epsilon=args.epsilon, cfg=cfg) - return gates + + print(gates) diff --git a/pygridsynth/config.py b/pygridsynth/config.py index 41a3287..421e522 100644 --- a/pygridsynth/config.py +++ b/pygridsynth/config.py @@ -11,7 +11,7 @@ class GridsynthConfig: floop: int = 10 dtimeout: float | None = None ftimeout: float | None = None - verbose: bool = False + verbose: int = 0 measure_time: bool = False show_graph: bool = False up_to_phase: bool = False diff --git a/pygridsynth/unitary.py b/pygridsynth/domega_unitary.py similarity index 50% rename from pygridsynth/unitary.py rename to pygridsynth/domega_unitary.py index 0406551..c871632 100644 --- a/pygridsynth/unitary.py +++ b/pygridsynth/domega_unitary.py @@ -1,11 +1,14 @@ from __future__ import annotations +import string from functools import cached_property import mpmath -from .quantum_gate import HGate, QuantumCircuit, SGate, SXGate, TGate, WGate -from .ring import DOmega +from .mymath import RealNum, einsum, from_matrix_to_tensor, from_tensor_to_matrix +from .quantum_circuit import QuantumCircuit +from .quantum_gate import CxGate, HGate, SGate, SingleQubitGate, SXGate, TGate, WGate +from .ring import OMEGA, DOmega class DOmegaUnitary: @@ -167,3 +170,151 @@ def from_gates(cls, gates: str) -> DOmegaUnitary: else: raise ValueError return unitary.reduce_denomexp() + + +class DOmegaMatrix: + def __init__( + self, mat: list[list[DOmega]], wires: list[int], k: int = 0, phase: RealNum = 0 + ) -> None: + n = len(wires) + if (len(mat) != 2**n) or any(len(row) != 2**n for row in mat): + raise ValueError( + f"Matrix must be a {2**n}x{2**n} square matrix to match wires" + f"(got {len(mat)}x{len(mat[0]) if len(mat) > 0 else 0})" + ) + + self._mat: list[list[DOmega]] = mat + self._wires: list[int] = wires + self._k: int = k + self._phase: mpmath.mpf = mpmath.mpf(phase) % (2 * mpmath.mp.pi) + + @property + def mat(self) -> list[list[DOmega]]: + return self._mat + + @property + def wires(self) -> list[int]: + return self._wires + + @property + def k(self) -> int: + return self._k + + @property + def phase(self) -> mpmath.mpf: + return self._phase + + @phase.setter + def phase(self, phase: RealNum) -> None: + self._phase = mpmath.mpf(phase) % (2 * mpmath.mp.pi) + + @classmethod + def from_domega_unitary( + cls, unitary: DOmegaUnitary, wires: list[int], phase: RealNum = 0 + ) -> DOmegaMatrix: + return DOmegaMatrix(unitary.to_matrix, phase=phase, wires=wires) + + @classmethod + def from_single_qubit_gate(cls, g: SingleQubitGate) -> DOmegaMatrix: + return DOmegaMatrix.from_domega_unitary( + DOmegaUnitary.from_circuit(QuantumCircuit.from_list([g])), + wires=[g.target_qubit], + ) + + @classmethod + def from_single_qubit_circuit( + cls, circuit: QuantumCircuit, wires: list[int] + ) -> DOmegaMatrix: + return DOmegaMatrix.from_domega_unitary( + DOmegaUnitary.from_circuit(circuit), phase=circuit.phase, wires=wires + ) + + @classmethod + def from_w_gate(cls, g: WGate) -> DOmegaMatrix: + return DOmegaMatrix([[DOmega(OMEGA, 0)]], wires=[]) + + @classmethod + def from_cx_gate(cls, g: CxGate) -> DOmegaMatrix: + wires = [g.control_qubit, g.target_qubit] + mat = [[DOmega.from_int(0) for j in range(4)] for i in range(4)] + for i, j in [(0, 0), (1, 1), (2, 3), (3, 2)]: + mat[i][j] = DOmega.from_int(1) + return DOmegaMatrix(mat, wires) + + @cached_property + def to_matrix(self) -> list[list[DOmega]]: + return self._mat + + @cached_property + def to_complex_matrix(self) -> mpmath.matrix: + return mpmath.matrix( + [[x.to_complex for x in row] for row in self._mat] + ) * mpmath.exp(1.0j * self._phase) + + def __repr__(self) -> str: + return ( + f"DOmegaMatrix({repr(self._mat)}, {self._wires}, {self._k}, {self._phase})" + ) + + def __eq__(self, other: DOmegaMatrix | object) -> bool: + if isinstance(other, DOmegaMatrix): + return self._mat == other.mat and self._k == other.k + else: + return False + + def __matmul__(self, other: DOmegaMatrix) -> DOmegaMatrix: + if isinstance(other, DOmegaMatrix): + all_wires = list(set(self._wires + other.wires)) + n_total = len(all_wires) + + labels = list(string.ascii_lowercase) + if (len(self._wires) + len(other.wires)) * 2 > len(labels): + raise ValueError("Too many qubits for automatic einsum labeling") + + out_row_labels = [labels[i] for i in range(n_total)] + out_col_labels = [labels[n_total + i] for i in range(n_total)] + + self_row_labels = [out_row_labels[all_wires.index(w)] for w in self._wires] + self_col_labels = [out_col_labels[all_wires.index(w)] for w in self._wires] + + other_row_labels = [out_row_labels[all_wires.index(w)] for w in other.wires] + other_col_labels = [out_col_labels[all_wires.index(w)] for w in other.wires] + + shared_wires = set(self._wires) & set(other.wires) + for idx, w in enumerate(shared_wires): + i = self._wires.index(w) + j = other.wires.index(w) + other_row_labels[j] = labels[n_total * 2 + idx] + self_col_labels[i] = labels[n_total * 2 + idx] + + einsum_str = ( + "".join(self_row_labels + self_col_labels) + + "," + + "".join(other_row_labels + other_col_labels) + + "->" + + "".join(out_row_labels + out_col_labels) + ) + + A_tensor = from_matrix_to_tensor(self._mat, len(self._wires)) + B_tensor = from_matrix_to_tensor(other.mat, len(other.wires)) + C_tensor = from_tensor_to_matrix( + einsum(einsum_str, A_tensor, B_tensor), n_total + ) + + return DOmegaMatrix( + C_tensor, + k=self._k + other.k, + phase=self._phase + other.phase, + wires=all_wires, + ) + else: + return NotImplemented + + @classmethod + def identity(cls, wires: list[int]) -> DOmegaMatrix: + n = 2 ** len(wires) + mat = [ + [DOmega.from_int(1) if i == j else DOmega.from_int(0) for j in range(n)] + for i in range(n) + ] + return cls(mat, wires) diff --git a/pygridsynth/gridsynth.py b/pygridsynth/gridsynth.py index 0c6560e..1a25d62 100644 --- a/pygridsynth/gridsynth.py +++ b/pygridsynth/gridsynth.py @@ -5,15 +5,23 @@ from .config import GridsynthConfig from .diophantine import Result, diophantine_dyadic +from .domega_unitary import DOmegaUnitary from .grid_op import GridOp -from .mymath import MPFConvertible, RealNum, solve_quadratic, sqrt -from .quantum_gate import QuantumCircuit, Rz +from .mymath import ( + MPFConvertible, + RealNum, + convert_theta_and_epsilon, + dps_for_epsilon, + solve_quadratic, + sqrt, +) +from .quantum_circuit import QuantumCircuit +from .quantum_gate import Rz from .region import ConvexSet, Ellipse, Rectangle from .ring import DOmega, DRootTwo, ZOmega, ZRootTwo from .synthesis_of_cliffordT import decompose_domega_unitary from .tdgp import solve_TDGP from .to_upright import to_upright_ellipse_pair, to_upright_set_pair -from .unitary import DOmegaUnitary class EpsilonRegion(ConvexSet): @@ -110,7 +118,7 @@ def error( "Either `dps` or `epsilon` must be provided to determine precision." ) else: - dps = _dps_for_epsilon(epsilon) + dps = dps_for_epsilon(epsilon) with mpmath.workdps(dps): theta = mpmath.mpf(theta) phase = mpmath.mpf(phase) @@ -149,7 +157,7 @@ def get_synthesized_unitary( "Either `dps` or `epsilon` must be provided to determine precision." ) else: - dps = _dps_for_epsilon(epsilon) + dps = dps_for_epsilon(epsilon) with mpmath.workdps(dps): return DOmegaUnitary.from_gates(gates).to_complex_matrix @@ -226,8 +234,8 @@ def _gridsynth_exact( tdgp_sets = (epsilon_region, unit_disk, *transformed) if cfg.measure_time: - print(f"to_upright_set_pair: {time.time() - start} s") - if cfg.verbose: + print(f"time of to_upright_set_pair: {(time.time() - start) * 1000} ms") + if cfg.verbose >= 2: print("------------------") u_approx = None @@ -253,7 +261,7 @@ def _gridsynth_exact( print( "time of diophantine_dyadic: " f"{time_of_diophantine_dyadic * 1000} ms" ) - if cfg.verbose: + if cfg.verbose >= 2: print(f"{u_approx=}") print("------------------") return u_approx @@ -290,8 +298,8 @@ def _gridsynth_up_to_phase( tdgp_sets1 = (epsilon_region1, unit_disk1, *transformed1) if cfg.measure_time: - print(f"to_upright_set_pair: {time.time() - start} s") - if cfg.verbose: + print(f"time of to_upright_set_pair: {(time.time() - start) * 1000} ms") + if cfg.verbose >= 2: print("------------------") u_approx = None @@ -330,7 +338,7 @@ def _gridsynth_up_to_phase( print( "time of diophantine_dyadic: " f"{time_of_diophantine_dyadic * 1000} ms" ) - if cfg.verbose: + if cfg.verbose >= 2: print(f"{u_approx=}") print("------------------") return u_approx @@ -345,41 +353,13 @@ def gridsynth( if cfg is None: cfg = GridsynthConfig(**kwargs) elif kwargs: - warnings.warn( - "When 'cfg' is provided, 'kwargs' are ignored.", - stacklevel=2, - ) + warnings.warn("When 'cfg' is provided, 'kwargs' are ignored.", stacklevel=2) if cfg.dps is None: - cfg.dps = _dps_for_epsilon(epsilon) - - if isinstance(theta, float): - warnings.warn( - ( - f"pygridsynth is synthesizing the angle {theta}. " - "Please verify that this is the intended value. " - "Using float may introduce precision errors; " - "consider using mpmath.mpf for exact precision." - ), - UserWarning, - stacklevel=2, - ) - - if isinstance(epsilon, float): - warnings.warn( - ( - f"pygridsynth is using epsilon={epsilon} as the tolerance. " - "Please verify that this is the intended value. " - "Using float may introduce precision errors; " - "consider using mpmath.mpf for exact precision." - ), - UserWarning, - stacklevel=2, - ) + cfg.dps = dps_for_epsilon(epsilon) with mpmath.workdps(cfg.dps): - theta = mpmath.mpf(theta) - epsilon = mpmath.mpf(epsilon) + theta, epsilon = convert_theta_and_epsilon(theta, epsilon, dps=cfg.dps) if cfg.up_to_phase: return _gridsynth_up_to_phase(theta, epsilon, cfg=cfg) @@ -397,13 +377,10 @@ def gridsynth_circuit( if cfg is None: cfg = GridsynthConfig(**kwargs) elif kwargs: - warnings.warn( - "When 'cfg' is provided, 'kwargs' are ignored.", - stacklevel=2, - ) + warnings.warn("When 'cfg' is provided, 'kwargs' are ignored.", stacklevel=2) if cfg.dps is None: - cfg.dps = _dps_for_epsilon(epsilon) + cfg.dps = dps_for_epsilon(epsilon) with mpmath.workdps(cfg.dps): start_total = time.time() if cfg.measure_time else 0.0 @@ -419,7 +396,7 @@ def gridsynth_circuit( print( f"time of decompose_domega_unitary: {(time.time() - start) * 1000} ms" ) - print(f"total time: {(time.time() - start_total) * 1000} ms") + print(f"time of gridsynth_circuit: {(time.time() - start_total) * 1000} ms") return circuit @@ -433,13 +410,10 @@ def gridsynth_gates( if cfg is None: cfg = GridsynthConfig(**kwargs) elif kwargs: - warnings.warn( - "When 'cfg' is provided, 'kwargs' are ignored.", - stacklevel=2, - ) + warnings.warn("When 'cfg' is provided, 'kwargs' are ignored.", stacklevel=2) if cfg.dps is None: - cfg.dps = _dps_for_epsilon(epsilon) + cfg.dps = dps_for_epsilon(epsilon) with mpmath.workdps(cfg.dps): circuit = gridsynth_circuit( @@ -449,9 +423,3 @@ def gridsynth_gates( cfg=cfg, ) return circuit.to_simple_str() - - -def _dps_for_epsilon(epsilon: MPFConvertible) -> int: - e = mpmath.mpf(epsilon) - k = -mpmath.log10(e) - return int(15 + 2.5 * int(mpmath.ceil(k))) # used in newsynth diff --git a/pygridsynth/mixed_synthesis.py b/pygridsynth/mixed_synthesis.py new file mode 100644 index 0000000..c7c37ab --- /dev/null +++ b/pygridsynth/mixed_synthesis.py @@ -0,0 +1,381 @@ +""" +Mixed synthesis module. +""" + +from __future__ import annotations + +import os +from multiprocessing import Pool + +import cvxpy as cp +import mpmath +import numpy as np +import scipy + +from .mixed_synthesis_utils import ( + get_random_hermitian_operator, + unitary_to_choi, + unitary_to_gptm, +) +from .multi_qubit_unitary_approximation import approximate_multi_qubit_unitary +from .mymath import dps_for_epsilon +from .quantum_circuit import QuantumCircuit +from .quantum_gate import CxGate, HGate, SGate, SXGate, TGate, WGate # noqa: F401 + + +def _get_default_solver() -> str | None: + """ + Get the default solver with priority order. + + Returns: + Available solver name, or None if no solver is available. + """ + available_solvers = cp.installed_solvers() + + preferred_solvers = ["GUROBI", "SCS", "ECOS", "OSQP"] + + for solver_name in preferred_solvers: + if solver_name in available_solvers: + return solver_name + + if available_solvers: + return available_solvers[0] + + return None + + +def solve_LP( + A: np.ndarray, + b: np.ndarray, + eps: float, + scale: float = 1.0, + solver: str | None = None, +) -> np.ndarray | None: + """ + Solve a linear programming problem for optimal mixing probabilities. + + Args: + A: Constraint matrix. + b: Constraint vector. + eps: Error tolerance parameter. + scale: Scale factor. + solver: Solver to use (None for auto-selection). + + Returns: + Optimal solution vector, or None on failure. + """ + n, m = A.shape + x = cp.Variable(m) + r = cp.Variable(n) + + x.value = np.ones(m) / (eps * m) + objective = cp.Minimize(cp.sum(r)) + constraints = [ + cp.sum(x) * eps == 1, + x >= 0, + r >= 0, + A @ x - b / eps <= r / scale, + -(A @ x - b / eps) <= r / scale, + ] + prob = cp.Problem(objective, constraints) + + if solver is None: + solver = _get_default_solver() + if solver is None: + return None + + try: + prob.solve(solver=solver, verbose=False) + except Exception: + return None + if x.value is None: + return None + else: + x = x.value / np.sum(x.value) + return x + + +def approximate_unitary_task( + eu: list[list[complex]], num_qubits: int, eps: float, dps: int = -1 +) -> tuple[float, list[str], np.ndarray, np.ndarray, np.ndarray]: + """ + Compute approximation task for a perturbed unitary matrix. + + Args: + eu: Target unitary matrix perturbed by a Hermitian operator (as a list). + num_qubits: Number of qubits. + eps: Error tolerance parameter. + dps: Decimal precision (default: -1 for auto). + + Returns: + Tuple of (phase, gates, eu_np, eu_gptm, eu_choi). + - phase: Phase of the circuit. + - gates: List of gate strings in the circuit. + - eu_np: Approximated unitary matrix (numpy array). + - eu_gptm: GPTM representation of approximated unitary. + - eu_choi: Choi representation of approximated unitary. + """ + if dps == -1: + dps = dps_for_epsilon(eps) + with mpmath.workdps(dps): + circuit, U_approx = approximate_multi_qubit_unitary( + mpmath.matrix(eu), num_qubits, mpmath.mpf(eps), return_domega_matrix=False + ) + eu_np = np.array(U_approx.tolist(), dtype=complex) + eu_gptm = unitary_to_gptm(eu_np) + eu_choi = unitary_to_choi(eu_np) + return float(circuit.phase), [str(g) for g in circuit], eu_np, eu_gptm, eu_choi + + +def compute_optimal_mixing_probabilities( + u_gptm: np.ndarray, + num_qubits: int, + eu_gptm_list: list[np.ndarray], + eu_choi_list: list[np.ndarray], + eps: float, +) -> tuple[np.ndarray, np.ndarray] | None: + """ + Compute optimal mixing probabilities + from GPTM representation of perturbed unitaries. + + Args: + u_gptm: GPTM representation of target unitary. + num_qubits: Number of qubits. + eu_gptm_list: List of GPTM representations of perturbed unitaries. + eu_choi_list: List of Choi representations of perturbed unitaries. + eps: Error tolerance parameter. + + Returns: + Tuple of (probs_gptm, u_choi_opt), or None on failure. + - probs_gptm: Array of mixing probabilities. + - u_choi_opt: Optimal Choi matrix of the mixed unitary. + """ + A = np.array([eu_gptm.reshape(eu_gptm.size) for eu_gptm in eu_gptm_list]).T.real + b = u_gptm.reshape(2 ** (4 * num_qubits)).real + + probs_gptm = solve_LP(A, b, eps=eps, scale=1e-2 / eps) + if probs_gptm is None: + return None + else: + u_choi_opt = np.einsum("i,ijk->jk", probs_gptm, eu_choi_list) + + return probs_gptm, u_choi_opt + + +def process_unitary_approximation_parallel( + unitary: mpmath.matrix, + num_qubits: int, + eps: float, + M: int, + seed: int = 123, + dps: int = -1, +) -> tuple[ + np.ndarray, + np.ndarray, + list[QuantumCircuit], + list[np.ndarray], + list[np.ndarray], + list[np.ndarray], +]: + """ + Process unitary approximation in parallel. + + Args: + unitary: Target unitary matrix. + num_qubits: Number of qubits. + eps: Error tolerance parameter. + M: Number of Hermitian operators for perturbation. + seed: Random seed. + dps: Decimal precision (default: -1 for auto). + + Returns: + Tuple of (u_gptm, u_choi, circuit_list, eu_np_list, eu_gptm_list, eu_choi_list). + - u_gptm: GPTM representation of target unitary. + - u_choi: Choi representation of target unitary. + - circuit_list: List of QuantumCircuit objects for perturbed unitaries. + - eu_np_list: List of target unitary matrices perturbed by Hermitian operators. + - eu_gptm_list: List of GPTM representations of perturbed unitaries. + - eu_choi_list: List of Choi representations of perturbed unitaries. + """ + herm_list = get_random_hermitian_operator(2**num_qubits, M, seed=seed) + unitary_np = np.array(unitary.tolist(), dtype=complex) + + u_choi = unitary_to_choi(unitary_np) + u_gptm = unitary_to_gptm(unitary_np) + eu_list = [ + (unitary @ mpmath.matrix(scipy.linalg.expm(1j * _herm * eps))).tolist() + for _herm in herm_list + ] + + num_workers = os.cpu_count() + with Pool(processes=num_workers) as pool: + results = pool.starmap( + approximate_unitary_task, + [(eu, num_qubits, eps, dps) for eu in eu_list], + ) + circuit_list = [ + QuantumCircuit(phase=result[0], args=[eval(g) for g in result[1]]) + for result in results + ] + eu_np_list = [result[2] for result in results] + eu_gptm_list = [result[3] for result in results] + eu_choi_list = [result[4] for result in results] + return (u_gptm, u_choi, circuit_list, eu_np_list, eu_gptm_list, eu_choi_list) + + +def process_unitary_approximation_sequential( + unitary: mpmath.matrix, + num_qubits: int, + eps: float, + M: int, + seed: int = 123, + dps: int = -1, +) -> tuple[ + np.ndarray, + np.ndarray, + list[QuantumCircuit], + list[np.ndarray], + list[np.ndarray], + list[np.ndarray], +]: + """ + Process unitary approximation sequentially. + + Args: + unitary: Target unitary matrix. + num_qubits: Number of qubits. + eps: Error tolerance parameter. + M: Number of Hermitian operators for perturbation. + seed: Random seed. + dps: Decimal precision (default: -1 for auto). + + Returns: + Tuple of (u_gptm, u_choi, circuit_list, eu_np_list, eu_gptm_list, eu_choi_list). + - u_gptm: GPTM representation of target unitary. + - u_choi: Choi representation of target unitary. + - circuit_list: List of QuantumCircuit objects for perturbed unitaries. + - eu_np_list: List of target unitary matrices perturbed by Hermitian operators. + - eu_gptm_list: List of GPTM representations of perturbed unitaries. + - eu_choi_list: List of Choi representations of perturbed unitaries. + """ + herm_list = get_random_hermitian_operator(2**num_qubits, M, seed=seed) + unitary_np = np.array(unitary.tolist(), dtype=complex) + u_choi = unitary_to_choi(unitary_np) + u_gptm = unitary_to_gptm(unitary_np) + eu_list = [ + (unitary @ mpmath.matrix(scipy.linalg.expm(1j * _herm * eps))).tolist() + for _herm in herm_list + ] + + circuit_list = [] + eu_np_list = [] + eu_gptm_list = [] + eu_choi_list = [] + for eu in eu_list: + result = approximate_unitary_task(eu, num_qubits, eps, dps=dps) + circuit = QuantumCircuit(phase=result[0], args=[eval(g) for g in result[1]]) + circuit_list.append(circuit) + eu_np_list.append(result[2]) + eu_gptm_list.append(result[3]) + eu_choi_list.append(result[4]) + return (u_gptm, u_choi, circuit_list, eu_np_list, eu_gptm_list, eu_choi_list) + + +def mixed_synthesis_parallel( + unitary: mpmath.matrix | np.ndarray, + num_qubits: int, + eps: float, + M: int, + seed: int = 123, + dps: int = -1, +) -> ( + tuple[list[QuantumCircuit], list[np.ndarray], np.ndarray, np.ndarray, np.ndarray] + | None +): + """ + Compute mixed probabilities for mixed unitary synthesis (parallel version). + + Args: + unitary: Target unitary matrix. + num_qubits: Number of qubits. + eps: Error tolerance parameter. + M: Number of Hermitian operators for perturbation. + seed: Random seed. + dps: Decimal precision (default: -1 for auto). + + Returns: + Tuple of (circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt), + or None on failure. + - circuit_list: List of QuantumCircuits for perturbed unitaries. + - eu_np_list: List of target unitary matrices perturbed by Hermitian operators. + - probs_gptm: Array of mixing probabilities. + - u_choi: Choi representation of target unitary. + - u_choi_opt: Optimal mixed Choi matrix. + """ + if not isinstance(unitary, mpmath.matrix): + unitary = mpmath.matrix(unitary) + + u_gptm, u_choi, circuit_list, eu_np_list, eu_gptm_list, eu_choi_list = ( + process_unitary_approximation_parallel( + unitary, num_qubits, eps, M, seed=seed, dps=dps + ) + ) + result = compute_optimal_mixing_probabilities( + u_gptm, num_qubits, eu_gptm_list, eu_choi_list, eps + ) + if result is None: + return None + + probs_gptm, u_choi_opt = result + + return circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt + + +def mixed_synthesis_sequential( + unitary: mpmath.matrix | np.ndarray, + num_qubits: int, + eps: float, + M: int, + seed: int = 123, + dps: int = -1, +) -> ( + tuple[list[QuantumCircuit], list[np.ndarray], np.ndarray, np.ndarray, np.ndarray] + | None +): + """ + Compute mixed probabilities for mixed unitary synthesis (sequential version). + + Args: + unitary: Target unitary matrix. + num_qubits: Number of qubits. + eps: Error tolerance parameter. + M: Number of Hermitian operators for perturbation. + seed: Random seed. + dps: Decimal precision (default: -1 for auto). + + Returns: + Tuple of (circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt), + or None on failure. + - circuit_list: List of QuantumCircuits for perturbed unitaries. + - eu_np_list: List of target unitary matrices perturbed by Hermitian operators. + - probs_gptm: Array of mixing probabilities. + - u_choi: Choi representation of target unitary. + - u_choi_opt: Optimal mixed Choi matrix. + """ + if not isinstance(unitary, mpmath.matrix): + unitary = mpmath.matrix(unitary) + + u_gptm, u_choi, circuit_list, eu_np_list, eu_gptm_list, eu_choi_list = ( + process_unitary_approximation_sequential( + unitary, num_qubits, eps, M, seed=seed, dps=dps + ) + ) + result = compute_optimal_mixing_probabilities( + u_gptm, num_qubits, eu_gptm_list, eu_choi_list, eps + ) + if result is None: + return None + + probs_gptm, u_choi_opt = result + + return circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt diff --git a/pygridsynth/mixed_synthesis_parallel.py b/pygridsynth/mixed_synthesis_parallel.py new file mode 100644 index 0000000..ec8a054 --- /dev/null +++ b/pygridsynth/mixed_synthesis_parallel.py @@ -0,0 +1,155 @@ +""" +Parallel execution version of mixed synthesis. +""" + +from __future__ import annotations + +import os +import warnings +from multiprocessing import Pool +from typing import TYPE_CHECKING + +import cvxpy as cp +import numpy as np + +from .mixed_synthesis import ( + compute_optimal_mixing_probabilities, + process_unitary_approximation_parallel, +) +from .mymath import diamond_norm_error_from_choi, random_su + +if TYPE_CHECKING: + from .quantum_circuit import QuantumCircuit + +warnings.filterwarnings("ignore", category=UserWarning) + + +def my_task( + args: tuple[ + int, np.ndarray, np.ndarray, int, list[np.ndarray], list[np.ndarray], float + ], +) -> tuple[int, np.ndarray | None, float | None]: + """ + Compute optimal mixing probabilities and diamond norm error for a task. + + Args: + args: Tuple of + (idx, u_gptm, u_choi, num_qubits, eu_gptm_list, eu_choi_list, eps). + - idx: Task index. + - u_gptm: GPTM representation of target unitary. + - u_choi: Choi representation of target unitary. + - num_qubits: Number of qubits. + - eu_gptm_list: List of GPTM representations of perturbed unitaries. + - eu_choi_list: List of Choi representations of perturbed unitaries. + - eps: Error tolerance parameter. + + Returns: + Tuple of (idx, probs_gptm, error). + - idx: Task index. + - probs_gptm: Array of mixing probabilities, or None on failure. + - error: Diamond norm error, or None on failure. + """ + idx, u_gptm, u_choi, num_qubits, eu_gptm_list, eu_choi_list, eps = args + result = compute_optimal_mixing_probabilities( + u_gptm, num_qubits, eu_gptm_list, eu_choi_list, eps + ) + if result is None: + return (idx, None, None) + probs_gptm, u_choi_opt = result + error = diamond_norm_error_from_choi(u_choi, u_choi_opt, eps, mixed_synthesis=True) + return (idx, probs_gptm, error) + + +def main() -> list[ + tuple[ + int, + float, + list[QuantumCircuit], + list[np.ndarray], + np.ndarray | None, + float | None, + ] +]: + """ + Main processing function. + + Returns: + List of results as tuples of + (num_qubits, eps, circuits, eu_np_list, probs_gptm, error). + - num_qubits: Number of qubits. + - eps: Error tolerance parameter. + - circuits: List of QuantumCircuits for perturbed unitaries. + - eu_np_list: List of approximated unitary matrices (numpy arrays). + - probs_gptm: Array of mixing probabilities, or None on failure. + - error: Diamond norm error, or None on failure. + """ + # Parameter settings + eps_list = [1e-3, 1e-4, 1e-5, 1e-6] + num_qubits_list = [1, 2] + M_list = [16, 64] + num_trial = 2 + + # Initialize GUROBI solver + cp.Problem(cp.Minimize(cp.Variable(1)), []).solve(solver=cp.GUROBI, verbose=False) + + final_results: list[ + tuple[ + int, + float, + list[QuantumCircuit], + list[np.ndarray], + np.ndarray | None, + float | None, + ] + ] = [] + tasks = [] + idx = 0 + + for num_qubits, M in zip(num_qubits_list, M_list): + unitaries = [random_su(num_qubits) for _ in range(num_trial)] + + # Process all unitary approximations first + + for eps in eps_list: + for unitary in unitaries: + # Process unitary approximation in parallel + u_gptm, u_choi, circuits, eu_np_list, eu_gptm_list, eu_choi_list = ( + process_unitary_approximation_parallel( + unitary, num_qubits, eps, M, seed=123 + ) + ) + final_results.append( + (num_qubits, eps, circuits, eu_np_list, None, None) + ) + tasks.append( + (idx, u_gptm, u_choi, num_qubits, eu_gptm_list, eu_choi_list, eps) + ) + idx += 1 + + # Compute optimal mixing probabilities in parallel + num_workers = os.cpu_count() + with Pool(processes=num_workers) as pool: + results = pool.map(my_task, tasks) + for result in results: + idx, probs_gptm, error = result + num_qubits, eps, circuits, eu_np_list, _, _ = final_results[idx] + final_results[idx] = ( + num_qubits, + eps, + circuits, + eu_np_list, + probs_gptm, + error, + ) + + return final_results + + +if __name__ == "__main__": + results = main() + for result in results: + num_qubits, eps, circuits, eu_np_list, probs_gptm, error = result + print(f"num_qubits: {num_qubits}") + print(f"eps: {eps}") + print(f"error: {error:.4e}") + print("--------------------------------") diff --git a/pygridsynth/mixed_synthesis_sequential.py b/pygridsynth/mixed_synthesis_sequential.py new file mode 100644 index 0000000..ab2a2f7 --- /dev/null +++ b/pygridsynth/mixed_synthesis_sequential.py @@ -0,0 +1,110 @@ +""" +Sequential execution version of mixed synthesis. +""" + +from __future__ import annotations + +import warnings +from typing import TYPE_CHECKING + +import cvxpy as cp +import numpy as np + +from .mixed_synthesis import ( + compute_optimal_mixing_probabilities, + process_unitary_approximation_sequential, +) +from .mymath import diamond_norm_error_from_choi, random_su + +if TYPE_CHECKING: + from .quantum_circuit import QuantumCircuit + +warnings.filterwarnings("ignore", category=UserWarning) + + +def main() -> list[ + tuple[ + int, + float, + list[QuantumCircuit], + list[np.ndarray], + np.ndarray | None, + float | None, + ] +]: + """ + Main processing function. + + Returns: + List of results as tuples of + (num_qubits, eps, circuits, eu_np_list, probs_gptm, error). + - num_qubits: Number of qubits. + - eps: Error tolerance parameter. + - circuits: List of QuantumCircuits for perturbed unitaries. + - eu_np_list: List of approximated unitary matrices (numpy arrays). + - probs_gptm: Array of mixing probabilities, or None on failure. + - error: Diamond norm error, or None on failure. + """ + # Parameter settings + eps_list = [1e-3, 1e-4] + num_qubits_list = [1, 2] + M_list = [16, 64] + num_trial = 2 + + # Initialize GUROBI solver + if "GUROBI" in cp.installed_solvers(): + cp.Problem(cp.Minimize(cp.Variable(1)), []).solve( + solver=cp.GUROBI, verbose=False + ) + + final_results: list[ + tuple[ + int, + float, + list[QuantumCircuit], + list[np.ndarray], + np.ndarray | None, + float | None, + ] + ] = [] + + for num_qubits, M in zip(num_qubits_list, M_list): + unitaries = [random_su(num_qubits) for _ in range(num_trial)] + + # Process each combination of eps and unitary + for eps in eps_list: + for unitary in unitaries: + # Process unitary approximation sequentially + u_gptm, u_choi, circuits, eu_np_list, eu_gptm_list, eu_choi_list = ( + process_unitary_approximation_sequential( + unitary, num_qubits, eps, M, seed=123 + ) + ) + # Compute optimal mixing probabilities + result = compute_optimal_mixing_probabilities( + u_gptm, num_qubits, eu_gptm_list, eu_choi_list, eps + ) + if result is None: + final_results.append( + (num_qubits, eps, circuits, eu_np_list, None, None) + ) + else: + probs_gptm, u_choi_opt = result + error = diamond_norm_error_from_choi( + u_choi, u_choi_opt, eps, mixed_synthesis=True + ) + final_results.append( + (num_qubits, eps, circuits, eu_np_list, probs_gptm, error) + ) + + return final_results + + +if __name__ == "__main__": + results = main() + for result in results: + num_qubits, eps, circuits, eu_np_list, probs_gptm, error = result + print(f"num_qubits: {num_qubits}") + print(f"eps: {eps}") + print(f"error: {error:.4e}") + print("--------------------------------") diff --git a/pygridsynth/mixed_synthesis_utils.py b/pygridsynth/mixed_synthesis_utils.py new file mode 100644 index 0000000..cf045b5 --- /dev/null +++ b/pygridsynth/mixed_synthesis_utils.py @@ -0,0 +1,488 @@ +""" +Utility functions for mixed synthesis. +""" + +from __future__ import annotations + +import cvxpy as cp +import numpy as np +from numba import njit + + +@njit +def vector_norm(x: np.ndarray) -> float: + """Compute Euclidean norm of a 1D array.""" + s = 0.0 + for i in range(x.shape[0]): + s += x[i] * x[i] + return np.sqrt(s) + + +@njit +def repulsive_update_numba( + points: np.ndarray, step: float, epsilon: float +) -> np.ndarray: + """ + Update points based on repulsive forces between them. + For each point i, compute repulsive force from all other points: + force = Σ_{j≠i} (x_i - x_j) / (||x_i - x_j||^2 + epsilon) + Project to tangent space, move by step, and project back to sphere. + """ + M, d = points.shape + new_points = np.empty_like(points) + for i in range(M): + force = np.zeros(d) + for j in range(M): + if i == j: + continue + diff = points[i] - points[j] + dist_sq = 0.0 + for k in range(d): + dist_sq += diff[k] * diff[k] + for k in range(d): + force[k] += diff[k] / (dist_sq + epsilon) + # Project to tangent space: subtract (force・x_i) x_i from force + dot = 0.0 + for k in range(d): + dot += force[k] * points[i][k] + tangent_force = np.empty(d) + for k in range(d): + tangent_force[k] = force[k] - dot * points[i][k] + # Update point and compute new position + x_new = np.empty(d) + for k in range(d): + x_new[k] = points[i][k] + step * tangent_force[k] + norm_val = vector_norm(x_new) + for k in range(d): + new_points[i][k] = x_new[k] / norm_val + return new_points + + +def relax_points_numba( + points: np.ndarray, + iterations: int = 100, + step: float = 0.001, + epsilon: float = 1e-6, +) -> np.ndarray: + """ + Relax point cloud using Numba-compiled repulsive_update_numba. + + Args: + points: Point cloud to relax. + iterations: Number of relaxation iterations. + step: Step size for updates. + epsilon: Small constant for numerical stability. + + Returns: + Relaxed point cloud. + """ + for _ in range(iterations): + points = repulsive_update_numba(points, step, epsilon) + return points + + +def initial_points(d: int, M: int, seed: int | None = None) -> np.ndarray: + """ + Uniformly sample M points from the unit sphere S^(d-1) in d dimensions. + + Args: + d: Dimension of the space. + M: Number of points to sample. + seed: Random seed (optional). + + Returns: + Array of shape (M, d) containing the sampled points. + """ + points = np.empty((M, d)) + if seed is not None: + np.random.seed(seed) + for i in range(M): + x = np.random.randn(d) + norm_val = np.linalg.norm(x) + points[i] = x / norm_val + return points + + +def vector_to_upper_triangular(vector: np.ndarray, d: int) -> np.ndarray: + """ + Convert a vector to an upper triangular matrix. + + Args: + vector: Vector to convert. + d: Dimension of the square matrix. + + Returns: + Upper triangular matrix (excluding diagonal). + """ + matrix = np.zeros((d, d)) + indices = np.triu_indices(d, k=1) + matrix[indices] = vector + return matrix + + +def points_to_hermitian_matrix(points: np.ndarray) -> np.ndarray: + """ + Convert points to a Hermitian matrix. + + Args: + points: Points array. + + Returns: + Hermitian matrix. + """ + hilbert_dim = int(np.sqrt(len(points) + 1)) + assert hilbert_dim**2 - 1 == len(points) + + non_diag_real = points[: (hilbert_dim * (hilbert_dim - 1)) // 2] + non_diag_imag = points[ + (hilbert_dim * (hilbert_dim - 1)) // 2 : hilbert_dim * (hilbert_dim - 1) + ] + diag = list(points[hilbert_dim * (hilbert_dim - 1) :]) + [ + -sum(points[(hilbert_dim * (hilbert_dim - 1)) :]) + ] + + mat = vector_to_upper_triangular( + non_diag_real, hilbert_dim + ) + 1j * vector_to_upper_triangular(non_diag_imag, hilbert_dim) + mat = mat + mat.conj().T + mat += np.diag(diag) + + return mat + + +def operate_M_tensor(vector: np.ndarray, M: np.ndarray) -> None: + r""" + In-place Fast Walsh-Hadamard Transform-like operation M^\otimes n @ vec. + + Args: + vector: Vector to transform (modified in-place). + M: Transformation matrix. + """ + h = 1 + d = M.shape[0] + + while h < len(vector): + # perform FWHT + for i in range(0, len(vector), h * d): + for j in range(i, i + h): + if d == 2: + x = vector[j % len(vector)] + y = vector[(j + h) % len(vector)] + vector[j % len(vector)] = M[0, 0] * x + M[0, 1] * y + vector[(j + h) % len(vector)] = M[1, 0] * x + M[1, 1] * y + elif d == 4: + x = vector[j % len(vector)] + y = vector[(j + h) % len(vector)] + z = vector[(j + 2 * h) % len(vector)] + w = vector[(j + 3 * h) % len(vector)] + + vector[j % len(vector)] = ( + M[0, 0] * x + M[0, 1] * y + M[0, 2] * z + M[0, 3] * w + ) + vector[(j + h) % len(vector)] = ( + M[1, 0] * x + M[1, 1] * y + M[1, 2] * z + M[1, 3] * w + ) + vector[(j + 2 * h) % len(vector)] = ( + M[2, 0] * x + M[2, 1] * y + M[2, 2] * z + M[2, 3] * w + ) + vector[(j + 3 * h) % len(vector)] = ( + M[3, 0] * x + M[3, 1] * y + M[3, 2] * z + M[3, 3] * w + ) + + # normalize and increment + h *= d + + +def permutate_qubit_vector(vector: np.ndarray, new_order: list[int]) -> np.ndarray: + """ + Change qubit order (replacement for qulacs permutate_qubit). + + Args: + vector: Vector of dimension 2^(n_qubit). + new_order: New qubit order (e.g., [0, 2, 1, 3] for 2 qubits). + + Returns: + Vector with reordered qubits. + """ + n_qubits = len(new_order) + dim = 2**n_qubits + + tensor = vector.reshape([2] * n_qubits) + transposed = np.transpose(tensor, new_order) + return transposed.reshape(dim) + + +def pauli_vec_to_state(pauli_vec: np.ndarray) -> np.ndarray: + """ + Convert Pauli vector to density matrix (implementation without qulacs). + + Args: + pauli_vec: Pauli vector. + + Returns: + Density matrix. + """ + n_qubit = (pauli_vec.shape[0].bit_length() - 1) // 2 + M_inv = np.array([[1, 0, 0, 1], [0, 1, -1j, 0], [0, 1, 1j, 0], [1, 0, 0, -1]]) + new_order = [2 * i for i in range(n_qubit)] + [2 * i + 1 for i in range(n_qubit)] + + rho_vec_tmp = pauli_vec.copy() + operate_M_tensor(rho_vec_tmp, M_inv) + + # Change qubit order (replacement for qulacs permutate_qubit) + rho_vec = permutate_qubit_vector(rho_vec_tmp, new_order) + + # Reshape to density matrix + rho = rho_vec.reshape(2**n_qubit, 2**n_qubit) + return rho + + +def get_random_hermitian_operator( + hilbert_dim: int, + num_herm: int = 1, + iterations: int = 100, + step: float = 0.001, + seed: int | None = None, +) -> list[np.ndarray]: + """ + Generate random Hermitian operators. + + Args: + hilbert_dim: Hilbert space dimension. + num_herm: Number of Hermitian operators to generate. + iterations: Number of relaxation iterations. + step: Step size for relaxation. + seed: Random seed (optional). + + Returns: + List of Hermitian matrices. + """ + d = hilbert_dim**2 - 1 + # Generate initial point cloud + points = initial_points(d, num_herm, seed=seed) + + # Lightweight relaxation using Numba + points_relaxed = relax_points_numba(points, iterations=iterations, step=step) + + herm_list = [] + for _points in points_relaxed: + if hilbert_dim == 2 ** (hilbert_dim.bit_length() - 1): + tmp = pauli_vec_to_state(np.concatenate(([0j], _points))) + herm = tmp / np.sqrt(2) + else: + tmp = points_to_hermitian_matrix(_points) + herm = tmp / np.sqrt(hilbert_dim) + herm_list.append(herm) + return herm_list + + +def hermitian_generalized_pauli_operators(d: int) -> list[np.ndarray]: + r""" + Construct Hermitian basis {H_k} (k=0,...,d^2-1) for d-dimensional Hilbert space + from Weyl-Heisenberg generalized Pauli operators Q_{a,b} = X^a Z^b, a,b = 0,...,d-1. + + Note: Standard Q_{a,b} are not Hermitian for d>2, so we construct Hermitian + basis from pairs Q_{a,b} and Q_{a,b}^\dagger: + H^{(R)}_{a,b} = (Q_{a,b}+Q_{a,b}^\dagger)/√2 + H^{(I)}_{a,b} = -i (Q_{a,b}-Q_{a,b}^\dagger)/√2 + + For self-adjoint cases (2a≡0, 2b≡0 mod d), we use phase correction: + H = exp(-iπ a b/d) Q_{a,b} + + Args: + d: Hilbert space dimension. + + Returns: + List of d^2 Hermitian (d,d) matrices (numpy arrays). + """ + H_ops = [] + omega = np.exp(2j * np.pi / d) + + # Define shift operator X + X = np.zeros((d, d), dtype=complex) + for j in range(d): + X[j, (j + 1) % d] = 1 + + # Define phase operator Z + Z = np.diag([omega**j for j in range(d)]) + + # Set for duplicate checking + processed = set() + + for a in range(d): + for b in range(d): + # Key for (a,b) (considering adjoint (a',b') for order independence) + key = (a, b) + # Adjoint corresponds to exponents (-a mod d, -b mod d) + a_p, b_p = (-a) % d, (-b) % d + key_partner = (a_p, b_p) + + # Skip if already processed + if key in processed or key_partner in processed: + continue + + # Q_{a,b} = X^a Z^b + Q = np.linalg.matrix_power(X, a) @ np.linalg.matrix_power(Z, b) + Q_dag = Q.conj().T # Q^\dagger + + # Self-adjoint case (up to phase): 2a≡0 and 2b≡0 (mod d) + if ((2 * a) % d == 0) and ((2 * b) % d == 0): + # Q_dag = ω^{-ab} Q, + # so phase correction φ = ω^{-ab/2} makes it Hermitian + phase = np.exp(-1j * np.pi * a * b / d) + H = phase * Q + H_ops.append(H) + processed.add(key) + else: + # Create Hermitian combinations from pairs + # Note: Q_dag = ω^{-ab} Q_{a_p,b_p}, + # but we use Q and Q_dag directly here + H_R = (Q + Q_dag) / np.sqrt(2) + H_I = -1j * (Q - Q_dag) / np.sqrt(2) + H_ops.append(H_R) + H_ops.append(H_I) + processed.add(key) + processed.add(key_partner) + + # Number of operators should be d^2 + if len(H_ops) != d**2: + raise ValueError( + f"Generated {len(H_ops)} Hermitian operators, but should be d^2 = {d**2}." + ) + + return H_ops + + +def unitary_to_gptm(U: np.ndarray) -> np.ndarray: + """ + Compute GPTM (Generalized Pauli Transfer Matrix) for unitary channel Λ(ρ) = U ρ U† + using Hermitian generalized Pauli basis {H_k}. + + Definition: + M_{ij} = (1/d) Tr[ H_j U H_i U† ] + + Args: + U: Unitary matrix of dimension (d,d). + + Returns: + GPTM matrix of size (d^2, d^2). + """ + d = U.shape[0] + H_ops = hermitian_generalized_pauli_operators(d) + num_ops = len(H_ops) + M = np.zeros((num_ops, num_ops), dtype=complex) + U_dag = U.conj().T + + for i, H_i in enumerate(H_ops): + transformed = U @ H_i @ U_dag + for j, H_j in enumerate(H_ops): + M[j, i] = np.trace(H_j @ transformed) / d + return M + + +def unitary_to_choi(unitary: np.ndarray) -> np.ndarray: + """ + Convert unitary matrix to Choi matrix. + + Args: + unitary: Unitary matrix. + + Returns: + Choi matrix. + """ + dim = unitary.shape[0] + max_entangled_vec = np.identity(dim).reshape(dim * dim) + + large_unitary = np.kron(np.identity(dim), unitary) + vec = large_unitary @ max_entangled_vec + + return np.outer(vec, vec.conj()) + + +def choi_to_unitary(choi: np.ndarray, tol: float = 1e-12) -> np.ndarray: + """ + Convert Choi matrix to unitary matrix. + + Args: + choi: Choi matrix. + tol: Tolerance for rank check. + + Returns: + Unitary matrix. + + Raises: + AssertionError: If input does not correspond to a unitary. + """ + vals, vecs = np.linalg.eigh(choi) + assert ( + np.linalg.matrix_rank(choi, tol=tol) == 1 + ), "input does not correspond to unitary." + n_qubit = (len(vals).bit_length() - 1) // 2 + + u_ret = vecs[:, np.argmax(vals)].reshape(2**n_qubit, 2**n_qubit).T * np.sqrt( + np.max(vals) + ) + u_ret /= u_ret[0, 0] / (np.abs(u_ret[0, 0])) + return u_ret + + +def diamond_norm_choi( + choi1: np.ndarray, + choi2: np.ndarray | None = None, + scale: float = 1, + solver: str | None = None, +) -> float: + """ + Compute the diamond norm of the difference channel whose Choi matrix is J_delta. + + Args: + choi1: Choi matrix of the first channel, shape (d*d, d*d). + choi2: Choi matrix of the second channel, shape (d*d, d*d) (optional). + scale: Scaling factor. + solver: The CVXPY solver to use (None for default: cp.SCS). + + Returns: + The computed diamond norm. + """ + if choi2 is None: + dim = int(np.sqrt(choi1.shape[0])) + choi2 = unitary_to_choi(np.diag(np.ones(dim))) + + if solver is None: + solver = cp.SCS + + J_delta = (choi1 - choi2) * scale + + d = int(np.sqrt(J_delta.shape[0])) + n = d * d + # Variable for the lifted operator in the SDP: an n x n Hermitian matrix. + X = cp.Variable((n, n), hermitian=True) + # Scalar variable t that will bound the operator norm of the partial trace. + t = cp.Variable(nonneg=True) + + constraints = [] + # X must be positive semidefinite and also dominate J_delta. + constraints += [X >> 0, X - J_delta >> 0] + + # We now define the partial trace of X over the second subsystem. + # For a block-structured matrix X (with blocks of size d x d), + # the (i,j) entry of Y = Tr_B(X) is given by + # summing the (i*d+k, j*d+k) entries of X. + Y = cp.Variable((d, d), hermitian=True) + for i in range(d): + for j in range(d): + # Sum over the appropriate indices to get (Y)_{ij} + constraints.append( + Y[i, j] == cp.sum([X[i * d + k, j * d + k] for k in range(d)]) + ) + + # Impose the operator norm constraint: ||Y||_∞ ≤ t. + # This is equivalent to: t*I - Y >= 0 and t*I + Y >= 0. + constraints.append(t * np.eye(d) - Y >> 0) + constraints.append(t * np.eye(d) + Y >> 0) + + # Set up the SDP: minimize t subject to the constraints. + prob = cp.Problem(cp.Minimize(t), constraints) + prob.solve(solver=solver) + + return prob.value * 2 / scale diff --git a/pygridsynth/multi_qubit_unitary_approximation.py b/pygridsynth/multi_qubit_unitary_approximation.py new file mode 100644 index 0000000..90ef50a --- /dev/null +++ b/pygridsynth/multi_qubit_unitary_approximation.py @@ -0,0 +1,339 @@ +import time +import warnings +from typing import Literal, overload + +import mpmath + +from .config import GridsynthConfig +from .domega_unitary import DOmegaMatrix +from .gridsynth import gridsynth_circuit +from .mymath import ( + MPFConvertible, + all_close, + convert_theta_and_epsilon, + dps_for_epsilon, + kron, +) +from .quantum_circuit import QuantumCircuit +from .quantum_gate import CxGate, HGate, QuantumGate, RzGate +from .two_qubit_unitary_approximation import approximate_two_qubit_unitary +from .unitary_approximation import approximate_one_qubit_unitary + + +def _blockZXZ( + U: mpmath.matrix, num_qubits: int, verbose: int = 0 +) -> tuple[mpmath.matrix, mpmath.matrix, mpmath.matrix, mpmath.matrix]: + n = 2**num_qubits + n_half = 2 ** (num_qubits - 1) + I = mpmath.eye(n_half) + + # upper left block + X = U[0:n_half, 0:n_half] + # lower left block + U12 = U[0:n_half, n_half:n] + + # svd: M = W*diag(S)*V.H + # X = W11*diag(S11)*Vh11 + # X = Sx*Ux + W11, S11, Vh11 = mpmath.svd_c(X) + Sx = W11 @ mpmath.diag(S11) @ W11.H + Ux = W11 @ Vh11 + + # svd: U12 = W12*diag(S12)*Vh12 + # Y = Sy*Uy + W12, S12, Vh12 = mpmath.svd_c(U12) + Sy = W12 @ mpmath.diag(S12) @ W12.H + Uy = W12 @ Vh12 + + A = (Sx + 1j * Sy) @ Ux + C = -1j * Ux.H @ Uy + B = U[n_half:n, 0:n_half] + U[n_half:n, n_half:n] @ C.H + Z = 2 * A.H @ X - I + + if verbose >= 1: + IC = mpmath.zeros(n) + IC[:n_half, :n_half] = I + IC[n_half:, n_half:] = C + + AB = mpmath.zeros(n) + AB[:n_half, :n_half] = A + AB[n_half:n, n_half:n] = B + + IZ = mpmath.zeros(n) + IZ[:n_half, :n_half] = I + Z + IZ[n_half:, n_half:] = I + Z + IZ[n_half:, :n_half] = I - Z + IZ[:n_half, n_half:] = I - Z + + print( + "Block-ZXZ correct for matrix of shape: ", + U.rows, + U.cols, + " ? ", + all_close(0.5 * AB @ IZ @ IC, U), + ) + + return A, B, Z, C + + +def _demultiplex( + M1: mpmath.matrix, M2: mpmath.matrix, verbose: int = 0 +) -> tuple[mpmath.matrix, list[mpmath.mpf], mpmath.matrix]: + eigenvalues, V = mpmath.eig(M1 @ M2.H) + D_sqrt = [mpmath.sqrt(eigenvalue) for eigenvalue in eigenvalues] + W = mpmath.diag(D_sqrt) @ V.H @ M2 + if verbose >= 1: + print( + "Demultiplexing correct? ", + all_close(V @ mpmath.diag(D_sqrt) @ W, M1), + all_close(V @ mpmath.diag(D_sqrt).conjugate() @ W, M2), + ) + return V, D_sqrt, W + + +def _bitwise_inner_product(a: int, b: int) -> int: + i = a & b + # number of set bits in i + i = i - ((i >> 1) & 0x55555555) + i = (i & 0x33333333) + ((i >> 2) & 0x33333333) + i = (((i + (i >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24 + return i % 2 + + +# returns M^k = (-1)^(b_(i-1)*g_(i-1)), +# where * is bitwise inner product, g = binary gray code, b = binary code. +def _genMk(n: int) -> mpmath.matrix: + Mk = mpmath.matrix(n) + for i in range(0, n): + for j in range(0, n): + Mk[i, j] = (-1) ** (_bitwise_inner_product((i), ((j) ^ (j) >> 1))) + return Mk + + +# input D as (one dimensional) array +def _decompose_cz( + D: mpmath.matrix, + control_qubits: list[int], + target_qubit: int, + inverse: bool = False, +) -> QuantumCircuit: + n = len(D) + control_qubits.reverse() + D_arg = [-2 * mpmath.arg(d) for d in D] + ar = mpmath.lu_solve(_genMk(n), D_arg) + circuit = QuantumCircuit() + + if n != 2 ** len(control_qubits): + print( + "Warning: shape mismatch for controlled Z between" + "# control bits and length of D" + ) + + if inverse: + for i in range(n - 1, -1, -1): + idx = (i ^ (i >> 1) ^ (i + 1) ^ ((i + 1) >> 1)).bit_length() - 1 + if idx == len(control_qubits): + posc = control_qubits[-1] + else: + posc = control_qubits[idx] + circuit.append(CxGate(posc, target_qubit)) + circuit.append(RzGate(ar[i].real, target_qubit)) + + else: + for i in range(0, n): + idx = (i ^ (i >> 1) ^ (i + 1) ^ ((i + 1) >> 1)).bit_length() - 1 + if idx == len(control_qubits): + posc = control_qubits[-1] + else: + posc = control_qubits[idx] + circuit.append(RzGate(ar[i].real, target_qubit)) + circuit.append(CxGate(posc, target_qubit)) + + return circuit + + +def _decompose_recursively( + U: mpmath.matrix, wires: list[int], verbose: int = 0 +) -> QuantumCircuit: + num_qubits = len(wires) + if num_qubits == 1 or num_qubits == 2: + return QuantumCircuit.from_list([QuantumGate(U, wires)]) + + n_half = 2 ** (num_qubits - 1) + I = mpmath.eye(n_half) + circuit = QuantumCircuit() + + A1, A2, B, C = _blockZXZ(U, num_qubits, verbose=verbose) + V_a, D_a, W_a = _demultiplex(A1, A2, verbose=verbose) + V_c, D_c, W_c = _demultiplex(I, C, verbose=verbose) + + V_a_circ = _decompose_recursively(V_a, wires[1:], verbose=verbose) + circuit += V_a_circ + + D_a_circ = _decompose_cz(D_a, wires[1:], wires[0]) + circuit += D_a_circ[:-1] + circuit.append(HGate(wires[0])) + + new_middle_ul = W_a @ V_c + Sz_I = kron(mpmath.matrix([[1, 0], [0, -1]]), mpmath.eye(2 ** (num_qubits - 2))) + new_middle_lr = Sz_I @ W_a @ B @ V_c @ Sz_I + V_m, D_sqrt_m, W_m = _demultiplex(new_middle_ul, new_middle_lr, verbose=verbose) + + V_m_circ = _decompose_recursively(V_m, wires[1:], verbose=verbose) + circuit += V_m_circ + + D_sqrt_m_circ = _decompose_cz(D_sqrt_m, wires[1:], wires[0]) + circuit += D_sqrt_m_circ + + W_m_circ = _decompose_recursively(W_m, wires[1:], verbose=verbose) + circuit += W_m_circ + + circuit.append(HGate(wires[0])) + + D_c_circ = _decompose_cz(D_c, wires[1:], wires[0], inverse=True) + circuit += D_c_circ[1:] + + circ_W_c = _decompose_recursively(W_c, wires[1:], verbose=verbose) + circuit += circ_W_c + + return circuit + + +@overload +def approximate_multi_qubit_unitary( + U: mpmath.matrix, + num_qubits: int, + epsilon: MPFConvertible, + cfg: GridsynthConfig | None = None, + return_domega_matrix: Literal[True] = True, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, DOmegaMatrix]: ... + + +@overload +def approximate_multi_qubit_unitary( + U: mpmath.matrix, + num_qubits: int, + epsilon: MPFConvertible, + cfg: GridsynthConfig | None = None, + return_domega_matrix: Literal[False] = False, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, mpmath.matrix]: ... + + +def approximate_multi_qubit_unitary( + U: mpmath.matrix, + num_qubits: int, + epsilon: MPFConvertible, + cfg: GridsynthConfig | None = None, + return_domega_matrix: bool = False, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, DOmegaMatrix | mpmath.matrix]: + if cfg is None: + cfg = GridsynthConfig(**kwargs) + elif kwargs: + warnings.warn("When 'cfg' is provided, 'kwargs' are ignored.", stacklevel=2) + + if cfg.dps is None: + cfg.dps = dps_for_epsilon(epsilon) + cfg.up_to_phase = True + + if num_qubits == 1: + return approximate_one_qubit_unitary( + U, + epsilon, + wires=[0], + decompose_partially=False, + return_domega_matrix=return_domega_matrix, + scale_epsilon=scale_epsilon, + cfg=cfg, + ) + elif num_qubits == 2: + return approximate_two_qubit_unitary( + U, + epsilon, + wires=[0, 1], + decompose_partially=False, + return_domega_matrix=return_domega_matrix, + scale_epsilon=scale_epsilon, + cfg=cfg, + ) + + with mpmath.workdps(cfg.dps): + _, epsilon = convert_theta_and_epsilon("0", epsilon, dps=cfg.dps) + if scale_epsilon: + epsilon /= 18 * 4 ** (num_qubits - 2) - 3 * 2 ** (num_qubits - 1) + 3 + wires = list(range(num_qubits)) + + start = time.time() if cfg.measure_time else 0.0 + decomposed_circuit = _decompose_recursively(U, wires=wires, verbose=cfg.verbose) + if cfg.measure_time: + print("--------------------------------") + print(f"time of _decompose_recursively: {(time.time() - start) * 1000} ms") + print("--------------------------------") + diag = [1] * (2**2) + circuit = QuantumCircuit() + + U_approx = DOmegaMatrix.identity(wires=wires) + for i, c in enumerate(decomposed_circuit): + if isinstance(c, RzGate): + circuit_rz = gridsynth_circuit(c.theta, epsilon, wires=c.wires, cfg=cfg) + circuit += circuit_rz + U_approx = U_approx @ DOmegaMatrix.from_single_qubit_circuit( + circuit_rz, wires=c.wires + ) + elif isinstance(c, CxGate): + circuit.append(c) + U_approx = U_approx @ DOmegaMatrix.from_cx_gate(c) + elif isinstance(c, HGate): + circuit.append(c) + U_approx = U_approx @ DOmegaMatrix.from_single_qubit_gate(c) + elif i == len(decomposed_circuit) - 1: + c.matrix = mpmath.diag(diag) @ c.matrix + scale = mpmath.det(c.matrix) ** (-1 / 4) + circuit.phase -= mpmath.arg(scale) + U_approx.phase -= mpmath.arg(scale) + c.matrix *= scale + tmp_circuit, tmp_unitary = approximate_two_qubit_unitary( + c.matrix, + epsilon, + wires=c.wires, + decompose_partially=False, + return_domega_matrix=True, + cfg=cfg, + ) + circuit += tmp_circuit + U_approx = U_approx @ tmp_unitary + else: + c.matrix = mpmath.diag(diag) @ c.matrix + scale = mpmath.det(c.matrix) ** (-1 / 4) + circuit.phase -= mpmath.arg(scale) + U_approx.phase -= mpmath.arg(scale) + c.matrix *= scale + tmp_circuit, tmp_diag, tmp_unitary = approximate_two_qubit_unitary( + c.matrix, + epsilon, + wires=c.wires, + decompose_partially=True, + return_domega_matrix=True, + scale_epsilon=False, + cfg=cfg, + ) + U_approx = U_approx @ tmp_unitary + circuit += tmp_circuit + diag = tmp_diag + + if cfg.measure_time: + print("--------------------------------") + print( + "time of approximate_multi_qubit_unitary: " + f"{(time.time() - start) * 1000} ms" + ) + print("--------------------------------") + if return_domega_matrix: + return circuit, U_approx + else: + return circuit, U_approx.to_complex_matrix diff --git a/pygridsynth/mymath.py b/pygridsynth/mymath.py index 0e73235..021656f 100644 --- a/pygridsynth/mymath.py +++ b/pygridsynth/mymath.py @@ -1,12 +1,56 @@ +import warnings from itertools import accumulate from typing import TypeAlias import mpmath +import numpy as np + +from .mixed_synthesis_utils import diamond_norm_choi, unitary_to_choi RealNum: TypeAlias = int | float | mpmath.mpf MPFConvertible: TypeAlias = RealNum | mpmath.mpf +def dps_for_epsilon(epsilon: MPFConvertible) -> int: + e = mpmath.mpf(epsilon) + k = -mpmath.log10(e) + return int(15 + 2.5 * int(mpmath.ceil(k))) # used in newsynth + + +def convert_theta_and_epsilon( + theta: MPFConvertible, epsilon: MPFConvertible, dps: int +) -> tuple[mpmath.mpf, mpmath.mpf]: + if isinstance(theta, float): + warnings.warn( + ( + f"pygridsynth is synthesizing the angle {theta}. " + "Please verify that this is the intended value. " + "Using float may introduce precision errors; " + "consider using mpmath.mpf for exact precision." + ), + UserWarning, + stacklevel=2, + ) + + if isinstance(epsilon, float): + warnings.warn( + ( + f"pygridsynth is using epsilon={epsilon} as the tolerance. " + "Please verify that this is the intended value. " + "Using float may introduce precision errors; " + "consider using mpmath.mpf for exact precision." + ), + UserWarning, + stacklevel=2, + ) + + with mpmath.workdps(dps): + theta_mpf = mpmath.mpf(theta) + epsilon_mpf = mpmath.mpf(epsilon) + + return theta_mpf, epsilon_mpf + + def SQRT2() -> mpmath.mpf: return mpmath.sqrt(2) @@ -99,3 +143,105 @@ def solve_quadratic( return (0, -b / a) else: return ((2 * c) / s2, (2 * c) / s1) + + +def mpmath_matrix_to_numpy(M: mpmath.matrix) -> np.ndarray: + return np.array(M.tolist(), dtype=complex) + + +def numpy_matrix_to_mpmath(M: np.ndarray) -> mpmath.matrix: + return mpmath.matrix(M.tolist()) + + +def trace(M: mpmath.matrix) -> mpmath.mpf: + return sum(M[i, i] for i in range(min(M.rows, M.cols))) + + +def kron(A: mpmath.matrix, B: mpmath.matrix) -> mpmath.matrix: + return mpmath.matrix(np.kron(A.tolist(), B.tolist())) + + +def all_close( + A: mpmath.matrix, B: mpmath.matrix, tol: mpmath.mpf = mpmath.mpf("1e-5") +) -> bool: + return mpmath.norm(mpmath.matrix(A - B), p="inf") < tol + + +def einsum(subscripts: str, *operands: np.ndarray) -> np.ndarray: + return np.einsum(subscripts, *operands) + + +def from_matrix_to_tensor(mat: list[list], n) -> np.ndarray: + return np.array(mat, dtype=object).reshape([2] * n * 2) + + +def from_tensor_to_matrix(mat: np.ndarray, n) -> list[list]: + return mat.reshape((2**n, 2**n)).tolist() + + +def random_su(n: int) -> mpmath.matrix: + """ + Generate a random SU(n) unitary matrix. + + Args: + n: Number of qubits. + + Returns: + Random SU(2^n) unitary matrix. + """ + dim = 2**n + z = mpmath.matrix(np.random.random_sample((dim, dim))) + 1j * mpmath.matrix( + np.random.random_sample((dim, dim)) + ) / mpmath.sqrt(2) + q, _ = mpmath.qr(z) + q /= mpmath.det(q) ** (1 / dim) + return q + + +def diamond_norm_error( + u: np.ndarray | mpmath.matrix, + u_opt: np.ndarray | mpmath.matrix, + eps: MPFConvertible, +) -> float: + """ + Compute error using diamond norm. + + Args: + u: Target unitary. + u_opt: Mixed unitary. + eps: Error tolerance parameter. + + Returns: + Diamond norm error between the target and mixed unitaries. + """ + if isinstance(u, mpmath.matrix): + u = mpmath_matrix_to_numpy(u) + if isinstance(u_opt, mpmath.matrix): + u_opt = mpmath_matrix_to_numpy(u_opt) + if isinstance(eps, mpmath.mpf): + eps = float(eps) + u_choi = unitary_to_choi(u) + u_choi_opt = unitary_to_choi(u_opt) + return diamond_norm_error_from_choi(u_choi, u_choi_opt, eps, mixed_synthesis=False) + + +def diamond_norm_error_from_choi( + u_choi: np.ndarray, + u_choi_opt: np.ndarray, + eps: float, + mixed_synthesis: bool = False, +) -> float: + """ + Compute error using diamond norm. + + Args: + u_choi: Choi representation of target unitary. + u_choi_opt: Choi representation of mixed unitary. + eps: Error tolerance parameter. + mixed_synthesis: Whether the error is for mixed synthesis. + + Returns: + Diamond norm error between the target and mixed unitaries. + """ + scale = 1e-2 / eps**2 if mixed_synthesis else 1e-2 / eps + return diamond_norm_choi(u_choi, u_choi_opt, scale=scale) diff --git a/pygridsynth/normal_form.py b/pygridsynth/normal_form.py index f6600e6..7389d29 100644 --- a/pygridsynth/normal_form.py +++ b/pygridsynth/normal_form.py @@ -5,15 +5,8 @@ import mpmath from .mymath import RealNum -from .quantum_gate import ( - HGate, - QuantumCircuit, - QuantumGate, - SGate, - SXGate, - TGate, - WGate, -) +from .quantum_circuit import QuantumCircuit +from .quantum_gate import HGate, QuantumGate, SGate, SXGate, TGate, WGate class Axis(Enum): diff --git a/pygridsynth/quantum_circuit.py b/pygridsynth/quantum_circuit.py new file mode 100644 index 0000000..f051677 --- /dev/null +++ b/pygridsynth/quantum_circuit.py @@ -0,0 +1,108 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from .domega_unitary import DOmegaMatrix + +import functools +import operator +from typing import Iterable + +import mpmath + +from .mymath import RealNum +from .quantum_gate import CxGate, QuantumGate, WGate, w_phase + + +class QuantumCircuit(list): + def __init__(self, phase: RealNum = 0, args: list[QuantumGate] = []) -> None: + self._phase: mpmath.mpf = mpmath.mpf(phase) % (2 * mpmath.mp.pi) + super().__init__(args) + + @classmethod + def from_list(cls, gates: list[QuantumGate]) -> QuantumCircuit: + return cls(args=gates) + + @property + def phase(self) -> mpmath.mpf: + return self._phase + + @phase.setter + def phase(self, phase: RealNum) -> None: + self._phase = mpmath.mpf(phase) % (2 * mpmath.mp.pi) + + def __str__(self) -> str: + return f"exp(1.j * {self._phase}) * " + " * ".join(str(gate) for gate in self) + + def to_simple_str(self) -> str: + return "".join(g.to_simple_str() for g in self) + + def to_domega_matrix(self, wires: list[int]) -> DOmegaMatrix: + from .domega_unitary import DOmegaMatrix + + product = [DOmegaMatrix.identity(wires=wires)] + product[-1].phase = self.phase + crt_matrix = None + for g in self: + if isinstance(g, WGate): + if crt_matrix is not None: + product.append(crt_matrix) + crt_matrix = None + product.append(DOmegaMatrix.from_w_gate(g)) + elif isinstance(g, CxGate): + if crt_matrix is not None: + product.append(crt_matrix) + crt_matrix = None + product.append(DOmegaMatrix.from_cx_gate(g)) + else: + if crt_matrix is None: + crt_matrix = DOmegaMatrix.from_single_qubit_gate(g) + else: + if crt_matrix.wires == g.wires: + crt_matrix = crt_matrix @ DOmegaMatrix.from_single_qubit_gate(g) + else: + product.append(crt_matrix) + crt_matrix = None + crt_matrix = DOmegaMatrix.from_single_qubit_gate(g) + if crt_matrix is not None: + product.append(crt_matrix) + return functools.reduce(operator.matmul, product) + + def to_complex_matrix(self, n: int) -> mpmath.matrix: + return self.to_domega_matrix(wires=list(range(n))).to_complex_matrix + + def __add__(self, other: QuantumCircuit | list) -> QuantumCircuit: + if isinstance(other, QuantumCircuit): + return QuantumCircuit(self.phase + other.phase, list(self) + list(other)) + elif isinstance(other, list): + return QuantumCircuit(self.phase, list(self) + other) + else: + return NotImplemented + + def __radd__(self, other: QuantumCircuit | list) -> QuantumCircuit: + if isinstance(other, QuantumCircuit): + return QuantumCircuit(other.phase + self.phase, list(other) + list(self)) + elif isinstance(other, list): + return QuantumCircuit(self.phase, other + list(self)) + else: + return NotImplemented + + def __iadd__(self, other: QuantumCircuit | Iterable) -> QuantumCircuit: + if isinstance(other, QuantumCircuit): + self.phase += other.phase + super().__iadd__(list(other)) + return self + elif isinstance(other, Iterable): + super().__iadd__(other) + return self + else: + return NotImplemented + + def decompose_phase_gate(self) -> None: + self._phase %= 2 * mpmath.mp.pi + + for _ in range(round(float(self._phase / w_phase()))): + self.append(WGate()) + + self._phase = 0 diff --git a/pygridsynth/quantum_gate.py b/pygridsynth/quantum_gate.py index 0ed22a0..736f4a6 100644 --- a/pygridsynth/quantum_gate.py +++ b/pygridsynth/quantum_gate.py @@ -1,16 +1,16 @@ from __future__ import annotations -from typing import Iterable - import mpmath -from .mymath import RealNum - def cnot01() -> mpmath.matrix: return mpmath.matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) +def cnot10() -> mpmath.matrix: + return mpmath.matrix([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]]) + + def w_phase() -> mpmath.mpf: return mpmath.mp.pi / 4 @@ -50,30 +50,6 @@ def wires(self) -> list[int]: def __str__(self) -> str: return f"QuantumGate({self.matrix}, {self.wires})" - def to_whole_matrix(self, num_qubits: int) -> mpmath.matrix: - whole_matrix = mpmath.matrix(2**num_qubits) - m = len(self._wires) - target_qubit_sorted = sorted(self._wires, reverse=True) - - for i in range(2**m): - for j in range(2**m): - k = 0 - while k < 2**num_qubits: - i2 = 0 - j2 = 0 - for l in range(m): - t = num_qubits - target_qubit_sorted[l] - 1 - if k & 1 << t: - k += 1 << t - i2 |= (i >> l & 1) << t - j2 |= (j >> l & 1) << t - if k >= 2**num_qubits: - break - whole_matrix[k | i2, k | j2] = self._matrix[i, j] - k += 1 - - return whole_matrix - class SingleQubitGate(QuantumGate): def __init__(self, matrix: mpmath.matrix, target_qubit: int) -> None: @@ -87,23 +63,6 @@ def target_qubit(self) -> int: def __str__(self) -> str: return f"SingleQubitGate({self.matrix}, {self.target_qubit})" - def to_whole_matrix(self, num_qubits: int) -> mpmath.matrix: - whole_matrix = mpmath.matrix(2**num_qubits) - t = num_qubits - self.target_qubit - 1 - - for i in range(2): - for j in range(2): - k = 0 - while k < 2**num_qubits: - if k & 1 << t: - k += 1 << t - if k >= 2**num_qubits: - break - whole_matrix[k | (i << t), k | (j << t)] = self._matrix[i, j] - k += 1 - - return whole_matrix - class HGate(SingleQubitGate): def __init__(self, target_qubit: int) -> None: @@ -152,9 +111,6 @@ def __str__(self) -> str: def to_simple_str(self) -> str: return "W" - def to_whole_matrix(self, num_qubits: int) -> mpmath.matrix: - return self._matrix - class SXGate(SingleQubitGate): def __init__(self, target_qubit: int) -> None: @@ -170,7 +126,7 @@ def to_simple_str(self) -> str: class RzGate(SingleQubitGate): def __init__(self, theta: mpmath.mpf, target_qubit: int) -> None: - self._theta = theta + self._theta: mpmath.mpf = theta matrix = Rz(theta) super().__init__(matrix, target_qubit) @@ -190,7 +146,7 @@ def __mul__(self, other: RzGate) -> RzGate: class RxGate(SingleQubitGate): def __init__(self, theta: mpmath.mpf, target_qubit: int) -> None: - self._theta = theta + self._theta: mpmath.mpf = theta matrix = Rx(theta) super().__init__(matrix, target_qubit) @@ -223,120 +179,3 @@ def target_qubit(self) -> int: def __str__(self) -> str: return f"CxGate({self.control_qubit}, {self.target_qubit})" - - def to_whole_matrix(self, num_qubits: int) -> mpmath.matrix: - whole_matrix = mpmath.matrix(2**num_qubits) - c = num_qubits - self.control_qubit - 1 - t = num_qubits - self.target_qubit - 1 - for i0 in range(2): - for j0 in range(2): - for i1 in range(2): - for j1 in range(2): - k = 0 - while k < 2**num_qubits: - if k & (1 << min(c, t)): - k += 1 << min(c, t) - if k & (1 << max(c, t)): - k += 1 << max(c, t) - if k >= 2**num_qubits: - break - i2 = k | (i0 << c) | (i1 << t) - j2 = k | (j0 << c) | (j1 << t) - whole_matrix[i2, j2] = self._matrix[ - i0 * 2 + i1, j0 * 2 + j1 - ] - k += 1 - return whole_matrix - - -class QuantumCircuit(list): - def __init__(self, phase: RealNum = 0, args: list[QuantumGate] = []) -> None: - self._phase = mpmath.mpf(phase) % (2 * mpmath.mp.pi) - super().__init__(args) - - @classmethod - def from_list(cls, gates: list[QuantumGate]) -> QuantumCircuit: - return cls(args=gates) - - @property - def phase(self) -> mpmath.mpf: - return self._phase - - @phase.setter - def phase(self, phase: RealNum) -> None: - self._phase = mpmath.mpf(phase) % (2 * mpmath.mp.pi) - - def __str__(self) -> str: - return f"exp(1.j * {self._phase}) * " + " * ".join(str(gate) for gate in self) - - def to_simple_str(self) -> str: - return "".join(g.to_simple_str() for g in self) - - def __add__(self, other: QuantumCircuit | list) -> QuantumCircuit: - if isinstance(other, QuantumCircuit): - return QuantumCircuit(self.phase + other.phase, list(self) + list(other)) - elif isinstance(other, list): - return QuantumCircuit(self.phase, list(self) + other) - else: - return NotImplemented - - def __radd__(self, other: QuantumCircuit | list) -> QuantumCircuit: - if isinstance(other, QuantumCircuit): - return QuantumCircuit(other.phase + self.phase, list(other) + list(self)) - elif isinstance(other, list): - return QuantumCircuit(self.phase, other + list(self)) - else: - return NotImplemented - - def __iadd__(self, other: QuantumCircuit | Iterable) -> QuantumCircuit: - if isinstance(other, QuantumCircuit): - self.phase += other.phase - super().__iadd__(list(other)) - return self - elif isinstance(other, Iterable): - super().__iadd__(other) - return self - else: - return NotImplemented - - def decompose_phase_gate(self) -> None: - self._phase %= 2 * mpmath.mp.pi - - for _ in range(round(float(self._phase / w_phase()))): - self.append(WGate()) - - self._phase = 0 - - def to_matrix(self, num_qubits: int) -> mpmath.matrix: - U = mpmath.eye(2**num_qubits) * mpmath.exp(1.0j * self._phase) - for g in self: - U2 = mpmath.zeros(2**num_qubits) - if isinstance(g, WGate): - U2 = U * g.matrix - elif isinstance(g, CxGate): - for i in range(2**num_qubits): - for j in range(2**num_qubits): - for b0 in range(2): - for b1 in range(2): - t = num_qubits - g.target_qubit - 1 - c = num_qubits - g.control_qubit - 1 - j2 = j ^ b0 << c ^ b1 << t - U2[i, j2] += ( - U[i, j] - * g.matrix[ - (j >> c & 1) << 1 | (j >> t & 1), - (j2 >> c & 1) << 1 | (j2 >> t & 1), - ] - ) - elif isinstance(g, SingleQubitGate): - t = num_qubits - g.target_qubit - 1 - for i in range(2**num_qubits): - for j in range(2**num_qubits): - U2[i, j] += U[i, j] * g.matrix[j >> t & 1, j >> t & 1] - j2 = j ^ 1 << t - U2[i, j2] += U[i, j] * g.matrix[j >> t & 1, j2 >> t & 1] - else: - raise ValueError - U = U2 - - return U diff --git a/pygridsynth/synthesis_of_cliffordT.py b/pygridsynth/synthesis_of_cliffordT.py index 08e9f89..c5b7a27 100644 --- a/pygridsynth/synthesis_of_cliffordT.py +++ b/pygridsynth/synthesis_of_cliffordT.py @@ -1,16 +1,8 @@ +from .domega_unitary import DOmegaUnitary from .normal_form import NormalForm -from .quantum_gate import ( - HGate, - QuantumCircuit, - QuantumGate, - SGate, - SXGate, - TGate, - WGate, - w_phase, -) +from .quantum_circuit import QuantumCircuit +from .quantum_gate import HGate, QuantumGate, SGate, SXGate, TGate, WGate, w_phase from .ring import OMEGA_POWER -from .unitary import DOmegaUnitary BIT_SHIFT = [0, 0, 1, 0, 2, 0, 1, 3, 3, 3, 0, 2, 2, 1, 0, 0] BIT_COUNT = [0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4] @@ -91,11 +83,11 @@ def decompose_domega_unitary( for _ in range(m_S): circuit.append(SGate(target_qubit=wires[0])) unitary = unitary.mul_by_S_power_from_left(-m_S) - if not up_to_phase: + if up_to_phase: + circuit.phase = m_W * w_phase() + else: for _ in range(m_W): circuit.append(WGate()) - else: - circuit.phase = m_W * w_phase() assert unitary == DOmegaUnitary.identity(), "decomposition failed..." circuit = NormalForm.from_circuit(circuit).to_circuit(wires=wires) diff --git a/pygridsynth/tdgp.py b/pygridsynth/tdgp.py index d4f9893..eb955b0 100644 --- a/pygridsynth/tdgp.py +++ b/pygridsynth/tdgp.py @@ -17,7 +17,7 @@ def solve_TDGP( bboxA: Rectangle, bboxB: Rectangle, k: int, - verbose: bool = False, + verbose: int = 0, show_graph: bool = False, ) -> Iterator[DOmega]: opG_inv = opG.inv diff --git a/pygridsynth/to_upright.py b/pygridsynth/to_upright.py index 5b57920..1ac8aa7 100644 --- a/pygridsynth/to_upright.py +++ b/pygridsynth/to_upright.py @@ -24,11 +24,11 @@ def _shift_ellipse_pair(ellipse_pair: EllipsePair, n: int) -> EllipsePair: def _step_lemma( - ellipse_pair: EllipsePair, opG_l: GridOp, opG_r: GridOp, verbose: bool = False + ellipse_pair: EllipsePair, opG_l: GridOp, opG_r: GridOp, verbose: int = 0 ) -> tuple[EllipsePair, GridOp, GridOp, bool]: A = ellipse_pair.A B = ellipse_pair.B - if verbose: + if verbose >= 3: print("-----") print(f"skew: {ellipse_pair.skew}, bias: {ellipse_pair.bias}") print( @@ -40,19 +40,19 @@ def _step_lemma( ) print("-----") if B.b < 0: - if verbose: + if verbose >= 3: print("Z") OP_Z = GridOp(ZOmega(0, 0, 0, 1), ZOmega(0, -1, 0, 0)) return _reduction(ellipse_pair, opG_l, opG_r, OP_Z) elif A.bias * B.bias < 1: - if verbose: + if verbose >= 3: print("X") OP_X = GridOp(ZOmega(0, 1, 0, 0), ZOmega(0, 0, 0, 1)) return _reduction(ellipse_pair, opG_l, opG_r, OP_X) elif ellipse_pair.bias > 33.971 or ellipse_pair.bias < 0.029437: n = round(log(ellipse_pair.bias) / log(LAMBDA.to_real) / 8) OP_S = GridOp(ZOmega(-1, 0, 1, 1), ZOmega(1, -1, 1, 0)) - if verbose: + if verbose >= 3: print(f"S ({n=})") return _reduction(ellipse_pair, opG_l, opG_r, OP_S**n) elif ellipse_pair.skew <= 15: @@ -60,7 +60,7 @@ def _step_lemma( elif ellipse_pair.bias > 5.8285 or ellipse_pair.bias < 0.17157: n = round(log(ellipse_pair.bias) / log(LAMBDA.to_real) / 4) ellipse_pair = _shift_ellipse_pair(ellipse_pair, n) - if verbose: + if verbose >= 3: print(f"sigma ({n=})") if n >= 0: OP_SIGMA_L = GridOp(ZOmega(-1, 0, 1, 1), ZOmega(0, 1, 0, 0)) ** n @@ -70,36 +70,36 @@ def _step_lemma( OP_SIGMA_R = GridOp(ZOmega(0, 0, 0, 1), ZOmega(1, 1, 1, 0)) ** (-n) return ellipse_pair, opG_l * OP_SIGMA_L, OP_SIGMA_R * opG_r, False elif 0.24410 <= A.bias <= 4.0968 and 0.24410 <= B.bias <= 4.0968: - if verbose: + if verbose >= 3: print("R") OP_R = GridOp(ZOmega(0, 0, 1, 0), ZOmega(1, 0, 0, 0)) return _reduction(ellipse_pair, opG_l, opG_r, OP_R) elif A.b >= 0 and A.bias <= 1.6969: - if verbose: + if verbose >= 3: print("K") OP_K = GridOp(ZOmega(-1, -1, 0, 0), ZOmega(0, -1, 1, 0)) return _reduction(ellipse_pair, opG_l, opG_r, OP_K) elif A.b >= 0 and B.bias <= 1.6969: - if verbose: + if verbose >= 3: print("K_conj_sq2") OP_K_conj_sq2 = GridOp(ZOmega(1, -1, 0, 0), ZOmega(0, -1, -1, 0)) return _reduction(ellipse_pair, opG_l, opG_r, OP_K_conj_sq2) elif A.b >= 0: n = max(1, floorsqrt(min(A.bias, B.bias) / 4)) - if verbose: + if verbose >= 3: print(f"A ({n=})") OP_A_n = GridOp(ZOmega(0, 0, 0, 1), ZOmega(0, 1, 0, 2 * n)) return _reduction(ellipse_pair, opG_l, opG_r, OP_A_n) else: n = max(1, floorsqrt(min(A.bias, B.bias) / 2)) - if verbose: + if verbose >= 3: print(f"B ({n=})") OP_B_n = GridOp(ZOmega(0, 0, 0, 1), ZOmega(n, 1, -n, 0)) return _reduction(ellipse_pair, opG_l, opG_r, OP_B_n) def to_upright_ellipse_pair( - ellipseA: Ellipse, ellipseB: Ellipse, verbose: bool = False + ellipseA: Ellipse, ellipseB: Ellipse, verbose: int = 0 ) -> GridOp: ellipseA_normalized = ellipseA.normalize() ellipseB_normalized = ellipseB.normalize() @@ -120,7 +120,7 @@ def to_upright_set_pair( setB: ConvexSet, opG: GridOp | None = None, show_graph: bool = False, - verbose: bool = False, + verbose: int = 0, ) -> tuple[GridOp, Ellipse, Ellipse, Rectangle, Rectangle]: if opG is None: opG = to_upright_ellipse_pair(setA.ellipse, setB.ellipse, verbose=verbose) @@ -131,7 +131,7 @@ def to_upright_set_pair( bboxB = ellipseB_upright.bbox() upA = ellipseA_upright.area / bboxA.area upB = ellipseB_upright.area / bboxB.area - if verbose: + if verbose >= 2: print(f"{upA=}, {upB=}") if show_graph: plot_sol([], ellipseA_upright, ellipseB_upright, bboxA, bboxB) diff --git a/pygridsynth/two_qubit_unitary_approximation.py b/pygridsynth/two_qubit_unitary_approximation.py new file mode 100644 index 0000000..85257cd --- /dev/null +++ b/pygridsynth/two_qubit_unitary_approximation.py @@ -0,0 +1,465 @@ +import time +import warnings +from typing import Literal, overload + +import mpmath +import numpy as np + +from .config import GridsynthConfig +from .domega_unitary import DOmegaMatrix +from .gridsynth import gridsynth_circuit +from .mymath import ( + MPFConvertible, + all_close, + convert_theta_and_epsilon, + dps_for_epsilon, + kron, + sqrt, + trace, +) +from .quantum_circuit import QuantumCircuit +from .quantum_gate import CxGate, HGate, Rx, RxGate, Rz, RzGate, SingleQubitGate, cnot01 +from .synthesis_of_cliffordT import decompose_domega_unitary +from .unitary_approximation import ( + approximate_one_qubit_unitary, + euler_decompose, + magnitude_approximate, +) + + +def _SU2SU2_to_tensor_products(U: mpmath.matrix) -> tuple[mpmath.matrix, mpmath.matrix]: + U_reshaped = mpmath.matrix( + np.array(U).reshape((2, 2, 2, 2)).transpose(0, 2, 1, 3).reshape(4, 4) + ) + u, s, vh = mpmath.svd_c(U_reshaped) + A = sqrt(s[0]) * mpmath.matrix(np.array(u[:, 0].tolist()).reshape((2, 2))) + B = sqrt(s[0]) * mpmath.matrix(np.array(vh[0, :].tolist()).reshape((2, 2))) + angle = mpmath.arg(mpmath.det(A)) / 2 + A *= mpmath.exp(-1.0j * angle) + B *= mpmath.exp(1.0j * angle) + return A, B + + +def _extract_SU2SU2_prefactors( + U: mpmath.matrix, V: mpmath.matrix +) -> tuple[mpmath.matrix, mpmath.matrix, mpmath.matrix, mpmath.matrix]: + u = E.H @ U @ E + v = E.H @ V @ E + uuT = u @ u.T + vvT = v @ v.T + uuT_real_plus_imag = mpmath.matrix( + [ + [mpmath.re(uuT[i, j]) + mpmath.im(uuT[i, j]) for j in range(4)] + for i in range(4) + ] + ) + vvT_real_plus_imag = mpmath.matrix( + [ + [mpmath.re(vvT[i, j]) + mpmath.im(vvT[i, j]) for j in range(4)] + for i in range(4) + ] + ) + _, p = mpmath.eigsy(uuT_real_plus_imag) + _, q = mpmath.eigsy(vvT_real_plus_imag) + p = p @ mpmath.diag([1, 1, 1, mpmath.sign(mpmath.det(p))]) + q = q @ mpmath.diag([1, 1, 1, mpmath.sign(mpmath.det(q))]) + + G = p @ q.T + H = v.H @ G.T @ u + AB = E @ G @ E.H + A, B = _SU2SU2_to_tensor_products(AB) + CD = E @ H @ E.H + C, D = _SU2SU2_to_tensor_products(CD) + + return A, B, C, D + + +Sy_Sy = mpmath.matrix([[0, 0, 0, -1], [0, 0, 1, 0], [0, 1, 0, 0], [-1, 0, 0, 0]]) +I_Sz = mpmath.matrix([[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]]) +E = ( + 1 + / sqrt(2) + * mpmath.matrix( + [[1, 1.0j, 0, 0], [0, 0, 1.0j, 1], [0, 0, 1.0j, -1], [1, -1.0j, 0, 0]] + ) +) + + +def _gamma(u: mpmath.matrix) -> mpmath.matrix: + return u @ Sy_Sy @ u.T @ Sy_Sy + + +def _decompose_with_2_cnots( + U: mpmath.matrix, wires: list[int], verbose: int = 0 +) -> QuantumCircuit: + m = _gamma(U) + + eig_vals, _ = mpmath.eig(m) + angle_eig_vals = sorted([mpmath.arg(eigval) for eigval in eig_vals]) + r = angle_eig_vals[0] + s = angle_eig_vals[1] + t = angle_eig_vals[2] + diff_r_s = mpmath.fabs(mpmath.exp(1.0j * r) + mpmath.exp(1.0j * s)) + diff_r_t = mpmath.fabs(mpmath.exp(1.0j * r) + mpmath.exp(1.0j * t)) + if diff_r_s < diff_r_t: + s, t = t, s + r = mpmath.fabs(r) + s = mpmath.fabs(s) + theta = (r + s) / 2 + phi = (r - s) / 2 + + V = cnot01() @ kron(Rx(theta), Rz(phi)) @ cnot01() + A, B, C, D = _extract_SU2SU2_prefactors(U, V) + + if verbose >= 1: + print( + "_decompose_with_2_cnots correct?", + all_close(U, kron(A, B) @ V @ kron(C, D)), + ) + + gates = [ + SingleQubitGate(A, wires[0]), + SingleQubitGate(B, wires[1]), + CxGate(wires[0], wires[1]), + RxGate(theta, wires[0]), + RzGate(phi, wires[1]), + CxGate(wires[0], wires[1]), + SingleQubitGate(C, wires[0]), + SingleQubitGate(D, wires[1]), + ] + circuit = QuantumCircuit.from_list(gates) + return circuit + + +@overload +def _decompose_two_qubit_unitary( + U: mpmath.matrix, + wires: list[int] = [0, 1], + decompose_partially: Literal[True] = True, + verbose: int = 0, +) -> tuple[QuantumCircuit, list[mpmath.mpf]]: ... + + +@overload +def _decompose_two_qubit_unitary( + U: mpmath.matrix, + wires: list[int] = [0, 1], + decompose_partially: Literal[False] = False, + verbose: int = 0, +) -> QuantumCircuit: ... + + +def _decompose_two_qubit_unitary( + U: mpmath.matrix, + wires: list[int] = [0, 1], + decompose_partially: bool = False, + verbose: int = 0, +) -> tuple[QuantumCircuit, list[mpmath.mpf]] | QuantumCircuit: + if decompose_partially: + gamma_U = _gamma(U.T) + + psi = mpmath.atan2(mpmath.im(trace(gamma_U)), mpmath.re(trace(I_Sz @ gamma_U))) + Delta = mpmath.diag( + [ + mpmath.exp(-1.0j * psi / 2), + mpmath.exp(1.0j * psi / 2), + mpmath.exp(1.0j * psi / 2), + mpmath.exp(-1.0j * psi / 2), + ] + ) + Delta_inv_diag = [ + mpmath.exp(1.0j * psi / 2), + mpmath.exp(-1.0j * psi / 2), + mpmath.exp(-1.0j * psi / 2), + mpmath.exp(1.0j * psi / 2), + ] + + circuit = _decompose_with_2_cnots(U @ Delta, wires=wires) + + if verbose >= 1: + print( + "_decompose_two_qubit_unitary correct?", + all_close( + U, + mpmath.exp(1.0j * circuit.phase) + * U + @ Delta + @ mpmath.diag(Delta_inv_diag), + ), + ) + return circuit, Delta_inv_diag + else: + U2 = mpmath.exp(0.25j * mpmath.mp.pi) * U @ cnot01() + gamma_U = _gamma(U2.T) + + psi = mpmath.atan2(mpmath.im(trace(gamma_U)), mpmath.re(trace(I_Sz @ gamma_U))) + + Delta = cnot01() @ kron(mpmath.eye(2), Rz(psi)) @ cnot01() + circuit = _decompose_with_2_cnots(U2 @ Delta, wires=wires) + + circuit.phase += -0.25 * mpmath.mp.pi + circuit += [CxGate(wires[0], wires[1]), RzGate(-psi, wires[1])] + if verbose >= 1: + print( + "_decompose_two_qubit_unitary correct?", + all_close( + U, + mpmath.exp(1.0j * circuit.phase) + * kron(circuit[0].matrix, circuit[1].matrix) + @ cnot01() + @ kron(circuit[3].matrix, circuit[4].matrix) + @ cnot01() + @ kron(circuit[6].matrix, circuit[7].matrix) + @ cnot01() + @ kron(mpmath.eye(2), circuit[9].matrix), + ), + ) + return circuit + + +@overload +def approximate_two_qubit_unitary( + U: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0, 1], + cfg: GridsynthConfig | None = None, + decompose_partially: Literal[True] = True, + return_domega_matrix: Literal[True] = True, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, list[mpmath.mpf], DOmegaMatrix]: ... + + +@overload +def approximate_two_qubit_unitary( + U: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0, 1], + cfg: GridsynthConfig | None = None, + decompose_partially: Literal[True] = True, + return_domega_matrix: Literal[False] = False, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, list[mpmath.mpf], mpmath.matrix]: ... + + +@overload +def approximate_two_qubit_unitary( + U: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0, 1], + cfg: GridsynthConfig | None = None, + decompose_partially: Literal[False] = False, + return_domega_matrix: Literal[True] = True, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, DOmegaMatrix]: ... + + +@overload +def approximate_two_qubit_unitary( + U: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0, 1], + cfg: GridsynthConfig | None = None, + decompose_partially: Literal[False] = False, + return_domega_matrix: Literal[False] = False, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, mpmath.matrix]: ... + + +@overload +def approximate_two_qubit_unitary( + U: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0, 1], + cfg: GridsynthConfig | None = None, + decompose_partially: Literal[False] = False, + return_domega_matrix: bool = False, + scale_epsilon: bool = True, + **kwargs, +) -> tuple[QuantumCircuit, DOmegaMatrix | mpmath.matrix]: ... + + +def approximate_two_qubit_unitary( + U: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0, 1], + cfg: GridsynthConfig | None = None, + decompose_partially: bool = False, + return_domega_matrix: bool = False, + scale_epsilon: bool = True, + **kwargs, +) -> ( + tuple[QuantumCircuit, list[mpmath.mpf], DOmegaMatrix | mpmath.matrix] + | tuple[QuantumCircuit, DOmegaMatrix | mpmath.matrix] +): + if cfg is None: + cfg = GridsynthConfig(**kwargs) + elif kwargs: + warnings.warn("When 'cfg' is provided, 'kwargs' are ignored.", stacklevel=2) + + if cfg.dps is None: + cfg.dps = dps_for_epsilon(epsilon) + cfg.up_to_phase |= decompose_partially + + with mpmath.workdps(cfg.dps): + _, epsilon = convert_theta_and_epsilon("0", epsilon, dps=cfg.dps) + if scale_epsilon: + if decompose_partially: + epsilon /= 12 + else: + epsilon /= 15 + + start = time.time() if cfg.measure_time else 0.0 + if decompose_partially: + decomposed_circuit, diag = _decompose_two_qubit_unitary( + U, + wires=wires, + decompose_partially=decompose_partially, + verbose=cfg.verbose, + ) + else: + decomposed_circuit = _decompose_two_qubit_unitary( + U, + wires=wires, + decompose_partially=decompose_partially, + verbose=cfg.verbose, + ) + if cfg.measure_time: + print("--------------------------------") + print( + "time of _decompose_two_qubit_unitary: " + f"{(time.time() - start) * 1000} ms" + ) + print("--------------------------------") + circuit = QuantumCircuit(phase=decomposed_circuit.phase) + rx_theta_approx = magnitude_approximate( + decomposed_circuit[3].theta, epsilon, cfg=cfg + ) + phase_mag_rx_theta, phi1_mag_rx_theta, theta_mag_rx_theta, phi2_mag_rx_theta = ( + euler_decompose(rx_theta_approx.to_complex_matrix) + ) + circuit_rx_theta_approx = decompose_domega_unitary( + rx_theta_approx, wires=decomposed_circuit[3].wires, up_to_phase=True + ) + rz_phi_approx = magnitude_approximate( + decomposed_circuit[4].theta, epsilon, cfg=cfg + ) + phase_mag_rz_phi, phi1_mag_rz_phi, theta_mag_rz_phi, phi2_mag_rz_phi = ( + euler_decompose(rz_phi_approx.to_complex_matrix) + ) + circuit_rz_phi_approx = decompose_domega_unitary( + rz_phi_approx, wires=decomposed_circuit[4].wires, up_to_phase=True + ) + circuit_rz_phi_approx = ( + [HGate(target_qubit=decomposed_circuit[4].target_qubit)] + + circuit_rz_phi_approx + + [HGate(target_qubit=decomposed_circuit[4].target_qubit)] + ) + + circuit.phase -= phase_mag_rx_theta + circuit.phase -= phase_mag_rz_phi + decomposed_circuit[0].matrix = decomposed_circuit[0].matrix @ Rz( + -phi1_mag_rx_theta + ) + decomposed_circuit[1].matrix = decomposed_circuit[1].matrix @ Rx( + -phi1_mag_rz_phi + ) + decomposed_circuit[6].matrix = ( + Rz(-phi2_mag_rx_theta) @ decomposed_circuit[6].matrix + ) + decomposed_circuit[7].matrix = ( + Rx(-phi2_mag_rz_phi) @ decomposed_circuit[7].matrix + ) + + U_approx = mpmath.eye(2**2) * mpmath.exp(1.0j * circuit.phase) + U_approx = DOmegaMatrix.identity(wires=wires) + U_approx.phase = circuit.phase + for i in range(len(decomposed_circuit)): + if isinstance(decomposed_circuit[i], CxGate): + cx_gate = CxGate( + control_qubit=decomposed_circuit[i].control_qubit, + target_qubit=decomposed_circuit[i].target_qubit, + ) + circuit.append(cx_gate) + U_approx = U_approx @ DOmegaMatrix.from_cx_gate(cx_gate) + else: + if i == 3: + circuit += circuit_rx_theta_approx + U_approx = U_approx @ DOmegaMatrix.from_single_qubit_circuit( + circuit_rx_theta_approx, [decomposed_circuit[i].target_qubit] + ) + elif i == 4: + circuit += circuit_rz_phi_approx + U_approx = U_approx @ DOmegaMatrix.from_single_qubit_circuit( + circuit_rz_phi_approx, [decomposed_circuit[i].target_qubit] + ) + elif decompose_partially and i in [6, 7]: + tmp_circuit, tmp_rz, tmp_unitary = approximate_one_qubit_unitary( + decomposed_circuit[i].matrix, + epsilon, + wires=decomposed_circuit[i].wires, + decompose_partially=decompose_partially, + return_domega_matrix=True, + scale_epsilon=False, + cfg=cfg, + ) + circuit += tmp_circuit + U_approx = U_approx @ tmp_unitary + + if decomposed_circuit[i].target_qubit == wires[0]: + diag[0] *= mpmath.exp(-1.0j * tmp_rz.theta / 2) + diag[1] *= mpmath.exp(-1.0j * tmp_rz.theta / 2) + diag[2] *= mpmath.exp(1.0j * tmp_rz.theta / 2) + diag[3] *= mpmath.exp(1.0j * tmp_rz.theta / 2) + elif decomposed_circuit[i].target_qubit == wires[1]: + diag[0] *= mpmath.exp(-1.0j * tmp_rz.theta / 2) + diag[1] *= mpmath.exp(1.0j * tmp_rz.theta / 2) + diag[2] *= mpmath.exp(-1.0j * tmp_rz.theta / 2) + diag[3] *= mpmath.exp(1.0j * tmp_rz.theta / 2) + elif i == 9: + circuit_rz_psi = gridsynth_circuit( + decomposed_circuit[i].theta, + epsilon, + wires=decomposed_circuit[i].wires, + cfg=cfg, + ) + circuit += circuit_rz_psi + U_approx = U_approx @ DOmegaMatrix.from_single_qubit_circuit( + circuit_rz_psi, [decomposed_circuit[i].target_qubit] + ) + + else: + tmp_circuit, tmp_unitary = approximate_one_qubit_unitary( + decomposed_circuit[i].matrix, + epsilon, + wires=decomposed_circuit[i].wires, + decompose_partially=False, + return_domega_matrix=True, + scale_epsilon=False, + cfg=cfg, + ) + circuit += tmp_circuit + U_approx = U_approx @ tmp_unitary + + if not cfg.up_to_phase: + circuit.decompose_phase_gate() + if cfg.measure_time: + print("--------------------------------") + print( + "time of approximate_two_qubit_unitary: " + f"{(time.time() - start) * 1000} ms" + ) + print("--------------------------------") + if decompose_partially: + if return_domega_matrix: + return circuit, diag, U_approx + else: + return circuit, diag, U_approx.to_complex_matrix + else: + if return_domega_matrix: + return circuit, U_approx + else: + return circuit, U_approx.to_complex_matrix diff --git a/pygridsynth/unitary_approximation.py b/pygridsynth/unitary_approximation.py new file mode 100644 index 0000000..6578125 --- /dev/null +++ b/pygridsynth/unitary_approximation.py @@ -0,0 +1,315 @@ +import time +import warnings +from typing import Literal, overload + +import mpmath + +from .config import GridsynthConfig +from .diophantine import Result, diophantine_dyadic +from .domega_unitary import DOmegaMatrix, DOmegaUnitary +from .gridsynth import gridsynth_circuit +from .mymath import MPFConvertible, convert_theta_and_epsilon, dps_for_epsilon, sqrt +from .normal_form import NormalForm +from .odgp import solve_scaled_ODGP +from .quantum_circuit import QuantumCircuit +from .quantum_gate import RzGate +from .region import Interval +from .ring import DRootTwo +from .synthesis_of_cliffordT import decompose_domega_unitary + + +def euler_decompose( + unitary: mpmath.matrix, +) -> tuple[mpmath.mpf, mpmath.mpf, mpmath.mpf, mpmath.mpf]: + # unitary = e^{i phase} e^{- i phi1 / 2 Z} e^{- i theta / 2 X} e^{- i phi2 / 2 Z} + # 0 <= theta <= pi, - pi <= phase < pi, - pi <= phi1 < pi, - pi <= phi2 < pi + + det_arg = mpmath.arg(mpmath.det(unitary)) + psi1 = mpmath.arg(unitary[1, 1]) + psi2 = mpmath.arg(unitary[1, 0]) + + phase = det_arg / 2 + theta = 2 * mpmath.atan2(mpmath.fabs(unitary[1, 0]), mpmath.fabs(unitary[1, 1])) + phi1 = psi1 + psi2 - det_arg + mpmath.mp.pi / 2 + phi2 = psi1 - psi2 - mpmath.mp.pi / 2 + + if phi1 >= mpmath.mp.pi: + phi1 -= 2 * mpmath.mp.pi + phase += mpmath.mp.pi + if phi1 < -mpmath.mp.pi: + phi1 += 2 * mpmath.mp.pi + phase += mpmath.mp.pi + if phi2 >= mpmath.mp.pi: + phi2 -= 2 * mpmath.mp.pi + phase += mpmath.mp.pi + if phi2 <= -mpmath.mp.pi: + phi2 += 2 * mpmath.mp.pi + phase += mpmath.mp.pi + if phase >= mpmath.mp.pi: + phase -= 2 * mpmath.mp.pi + + return (phase, phi1, theta, phi2) + + +def _generate_epsilon_interval(theta: mpmath.mpf, epsilon: mpmath.mpf) -> Interval: + l = mpmath.cos(-theta / 2 - epsilon / 2) ** 2 + r = mpmath.cos(-theta / 2 + epsilon / 2) ** 2 + if l > r: + l, r = r, l + if -epsilon <= theta <= epsilon: + r = 1 + if mpmath.mp.pi - epsilon / 2 <= theta / 2 <= mpmath.mp.pi + epsilon / 2: + r = 1 + if -mpmath.mp.pi - epsilon / 2 <= theta / 2 <= -mpmath.mp.pi + epsilon / 2: + r = 1 + if mpmath.mp.pi / 2 - epsilon / 2 <= theta / 2 <= mpmath.mp.pi / 2 + epsilon / 2: + l = 0 + + return Interval(l, r) + + +def unitary_diamond_distance(u1: mpmath.matrix, u2: mpmath.matrix) -> mpmath.mpf: + mat = u1.H @ u2 + eig_vals, _ = mpmath.eig(mat) + phases = [mpmath.arg(eig_val) for eig_val in eig_vals] + d = mpmath.fabs(mpmath.cos((phases[0] - phases[1]) / 2)) + return 2 * sqrt(1 - d**2) + + +def check(u_target: mpmath.matrix, gates: str) -> None: + t_count = gates.count("T") + h_count = gates.count("H") + u_approx = DOmegaUnitary.from_gates(gates) + e = unitary_diamond_distance(u_target, u_approx.to_complex_matrix) + print(f"{gates=}") + print(f"{t_count=}, {h_count=}") + print(f"u_approx={u_approx.to_matrix}") + print(f"{e=}") + + +def _magnitude_approximate_with_fixed_k( + odgp_sets: tuple[Interval, Interval], + k: int, + cfg: GridsynthConfig, +) -> tuple[DOmegaUnitary | None, float, float]: + start = time.time() if cfg.measure_time else 0.0 + sol = solve_scaled_ODGP(*odgp_sets, k) + time_of_solve_ODGP = time.time() - start if cfg.measure_time else 0.0 + + start = time.time() if cfg.measure_time else 0.0 + time_of_diophantine_dyadic = 0.0 + u_approx = None + for m in sol: + if m.parity == 0: + continue + z = diophantine_dyadic(m, seed=cfg.seed, loop_controller=cfg.loop_controller) + if not isinstance(z, Result): + if (z * z.conj).residue == 0: + continue + xi = 1 - DRootTwo.fromDOmega(z.conj * z) + w = diophantine_dyadic( + xi, seed=cfg.seed, loop_controller=cfg.loop_controller + ) + if not isinstance(w, Result): + z = z.reduce_denomexp() + w = w.reduce_denomexp() + if z.k > w.k: + w = w.renew_denomexp(z.k) + elif z.k < w.k: + z = z.renew_denomexp(w.k) + + k1 = (z + w).reduce_denomexp().k + k2 = (z + w.mul_by_omega()).reduce_denomexp().k + if k1 <= k2: + u_approx = DOmegaUnitary(z, w, 0) + else: + u_approx = DOmegaUnitary(z, w.mul_by_omega(), 0) + break + + time_of_diophantine_dyadic += time.time() - start if cfg.measure_time else 0.0 + return u_approx, time_of_solve_ODGP, time_of_diophantine_dyadic + + +def magnitude_approximate( + theta: MPFConvertible, + epsilon: MPFConvertible, + cfg: GridsynthConfig | None = None, + **kwargs, +) -> DOmegaUnitary: + if cfg is None: + cfg = GridsynthConfig(**kwargs) + elif kwargs: + warnings.warn( + "When 'cfg' is provided, 'kwargs' are ignored.", + stacklevel=2, + ) + + if cfg.dps is None: + cfg.dps = dps_for_epsilon(epsilon) + + with mpmath.workdps(cfg.dps): + theta, epsilon = convert_theta_and_epsilon(theta, epsilon, dps=cfg.dps) + + epsilon_interval = _generate_epsilon_interval(theta, epsilon) + unit_interval = Interval(0, 1) + odgp_sets = (epsilon_interval, unit_interval) + + u_approx = None + k = 0 + time_of_solve_ODGP = 0.0 + time_of_diophantine_dyadic = 0.0 + while True: + u_approx, time1, time2 = _magnitude_approximate_with_fixed_k( + odgp_sets, k, cfg=cfg + ) + time_of_solve_ODGP += time1 + time_of_diophantine_dyadic += time2 + if u_approx is not None: + break + + k += 1 + + if cfg.measure_time: + print(f"time of solve_ODGP: {time_of_solve_ODGP * 1000} ms") + print( + "time of diophantine_dyadic: " f"{time_of_diophantine_dyadic * 1000} ms" + ) + if cfg.verbose >= 2: + print(f"{u_approx=}") + return u_approx + + +@overload +def approximate_one_qubit_unitary( + unitary: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0], + decompose_partially: Literal[True] = True, + return_domega_matrix: Literal[True] = True, + scale_epsilon: bool = True, + cfg: GridsynthConfig | None = None, + **kwargs, +) -> tuple[QuantumCircuit, RzGate, DOmegaMatrix]: ... + + +@overload +def approximate_one_qubit_unitary( + unitary: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0], + decompose_partially: Literal[True] = True, + return_domega_matrix: Literal[False] = False, + scale_epsilon: bool = True, + cfg: GridsynthConfig | None = None, + **kwargs, +) -> tuple[QuantumCircuit, RzGate, mpmath.matrix]: ... + + +@overload +def approximate_one_qubit_unitary( + unitary: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0], + decompose_partially: Literal[False] = False, + return_domega_matrix: Literal[True] = True, + scale_epsilon: bool = True, + cfg: GridsynthConfig | None = None, + **kwargs, +) -> tuple[QuantumCircuit, DOmegaMatrix]: ... + + +@overload +def approximate_one_qubit_unitary( + unitary: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0], + decompose_partially: Literal[False] = False, + return_domega_matrix: Literal[False] = False, + scale_epsilon: bool = True, + cfg: GridsynthConfig | None = None, + **kwargs, +) -> tuple[QuantumCircuit, mpmath.matrix]: ... + + +@overload +def approximate_one_qubit_unitary( + unitary: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0], + decompose_partially: Literal[False] = False, + return_domega_matrix: bool = False, + scale_epsilon: bool = True, + cfg: GridsynthConfig | None = None, + **kwargs, +) -> tuple[QuantumCircuit, DOmegaMatrix | mpmath.matrix]: ... + + +def approximate_one_qubit_unitary( + unitary: mpmath.matrix, + epsilon: MPFConvertible, + wires: list[int] = [0], + decompose_partially: bool = False, + return_domega_matrix: bool = False, + scale_epsilon: bool = True, + cfg: GridsynthConfig | None = None, + **kwargs, +) -> ( + tuple[QuantumCircuit, RzGate, DOmegaMatrix | mpmath.matrix] + | tuple[QuantumCircuit, DOmegaMatrix | mpmath.matrix] +): + if cfg is None: + cfg = GridsynthConfig(**kwargs) + elif kwargs: + warnings.warn("When 'cfg' is provided, 'kwargs' are ignored.", stacklevel=2) + + if cfg.dps is None: + cfg.dps = dps_for_epsilon(epsilon) + cfg.up_to_phase |= decompose_partially + + with mpmath.workdps(cfg.dps): + _, epsilon = convert_theta_and_epsilon("0", epsilon, dps=cfg.dps) + if scale_epsilon: + if decompose_partially: + epsilon /= 2 + else: + epsilon /= 3 + + phase, phi1, theta, phi2 = euler_decompose(unitary) + rx_approx = magnitude_approximate(theta, epsilon, cfg=cfg) + phase_mag, phi1_mag, theta_mag, phi2_mag = euler_decompose( + rx_approx.to_complex_matrix + ) + phase -= phase_mag + phi1 -= phi1_mag + phi2 -= phi2_mag + + circuit_rz_left = gridsynth_circuit(phi1, epsilon, wires=wires, cfg=cfg) + circuit_rx_approx = decompose_domega_unitary( + rx_approx, wires=wires, up_to_phase=cfg.up_to_phase + ) + circuit = circuit_rz_left + circuit_rx_approx + circuit.phase += phase + Rz_right = RzGate(phi2, wires[0]) + + if decompose_partially: + circuit = NormalForm.from_circuit(circuit).to_circuit(wires=wires) + U_approx = DOmegaMatrix.from_single_qubit_circuit(circuit, wires=wires) + if return_domega_matrix: + return circuit, Rz_right, U_approx + else: + return circuit, Rz_right, U_approx.to_complex_matrix + + else: + circuit_rz_right = gridsynth_circuit( + Rz_right.theta, epsilon, wires=wires, cfg=cfg + ) + circuit += circuit_rz_right + + if not cfg.up_to_phase: + circuit.decompose_phase_gate() + circuit = NormalForm.from_circuit(circuit).to_circuit(wires=wires) + U_approx = DOmegaMatrix.from_single_qubit_circuit(circuit, wires=wires) + if return_domega_matrix: + return circuit, U_approx + else: + return circuit, U_approx.to_complex_matrix diff --git a/pyproject.toml b/pyproject.toml index d239f20..9e716cc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,6 +16,10 @@ dependencies = [ "mpmath>=1.3", "setuptools", "matplotlib", + "numpy", + "scipy", + "cvxpy", + "numba", ] [project.optional-dependencies] @@ -50,7 +54,7 @@ profile = "black" [tool.flake8] max-line-length = 88 -extend-ignore = "E203, E741, W503" +extend-ignore = "E203, E704, E741, W503" [dependency-groups] dev = [ diff --git a/requirements-dev.txt b/requirements-dev.txt index 215ac55..2b8f551 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,2 +1,5 @@ -mpmath +mpmath>=1.3 numpy +scipy +cvxpy +numba diff --git a/tests/test_main.py b/tests/test_main.py index 0570615..8001c85 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -6,7 +6,7 @@ from pygridsynth.__main__ import main -def test_main_deterministic_results(): +def test_main_deterministic_results(capsys): """Test that main function returns the same result for identical inputs""" # pi/8 ≈ 0.39269908169 test_args = ["pygridsynth", "0.39269908169", "0.01"] @@ -16,8 +16,9 @@ def test_main_deterministic_results(): with patch("sys.argv", test_args): original_dps = mpmath.mp.dps try: - result = main() - results.append(str(result)) + main() + captured = capsys.readouterr() + results.append(captured.out.strip()) finally: mpmath.mp.dps = original_dps @@ -25,7 +26,7 @@ def test_main_deterministic_results(): print(f"✓ Got the same result from 5 executions: {results[0]}") -def test_main_with_different_inputs(): +def test_main_with_different_inputs(capsys): """Test that main function works correctly with different inputs""" test_cases = [ ["pygridsynth", "0.78539816339", "0.1"], # pi/4 @@ -38,7 +39,9 @@ def test_main_with_different_inputs(): with patch("sys.argv", test_args): original_dps = mpmath.mp.dps try: - result = main() + main() + captured = capsys.readouterr() + result = captured.out.strip() results.append(result) assert result is not None, f"Result is None for args {test_args}" finally: @@ -47,9 +50,17 @@ def test_main_with_different_inputs(): print(f"✓ Successfully executed with {len(test_cases)} different inputs") -def test_main_consistency_with_options(): +def test_main_consistency_with_options(capsys): """Test that consistent results are obtained even with options""" - test_args = ["pygridsynth", "0.39269908169", "0.01", "--dps", "30", "--verbose"] + test_args = [ + "pygridsynth", + "0.39269908169", + "0.01", + "--dps", + "30", + "--verbose", + "1", + ] results = [] # Execute 3 times with the same arguments @@ -57,8 +68,9 @@ def test_main_consistency_with_options(): with patch("sys.argv", test_args): original_dps = mpmath.mp.dps try: - result = main() - results.append(str(result)) + main() + captured = capsys.readouterr() + results.append(captured.out.strip()) finally: mpmath.mp.dps = original_dps diff --git a/tests/test_mixed_synthesis.py b/tests/test_mixed_synthesis.py new file mode 100644 index 0000000..208f81e --- /dev/null +++ b/tests/test_mixed_synthesis.py @@ -0,0 +1,190 @@ +import numpy as np +import pytest + +from pygridsynth.mixed_synthesis import ( + mixed_synthesis_parallel, + mixed_synthesis_sequential, +) +from pygridsynth.mymath import random_su + + +def test_mixed_synthesis_sequential_basic(): + """Test basic functionality of mixed_synthesis_sequential""" + num_qubits = 2 + eps = 1e-4 + M = 64 + unitary = random_su(num_qubits) + + result = mixed_synthesis_sequential(unitary, num_qubits, eps, M, seed=123) + + assert result is not None, "Result should not be None" + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + + assert isinstance(circuit_list, list), "circuit_list should be a list" + assert len(circuit_list) == M, f"circuit_list should have {M} elements" + assert isinstance(eu_np_list, list), "eu_np_list should be a list" + assert len(eu_np_list) == M, f"eu_np_list should have {M} elements" + assert isinstance(probs_gptm, np.ndarray), "probs_gptm should be a numpy array" + assert len(probs_gptm) == M, f"probs_gptm should have {M} elements" + assert isinstance(u_choi, np.ndarray), "u_choi should be a numpy array" + assert isinstance(u_choi_opt, np.ndarray), "u_choi_opt should be a numpy array" + assert np.allclose(np.sum(probs_gptm), 1.0), "Probabilities should sum to 1" + assert np.all(probs_gptm >= 0), "All probabilities should be non-negative" + + +def test_mixed_synthesis_parallel_basic(): + """Test basic functionality of mixed_synthesis_parallel""" + num_qubits = 2 + eps = 1e-4 + M = 64 + unitary = random_su(num_qubits) + + result = mixed_synthesis_parallel(unitary, num_qubits, eps, M, seed=123) + + assert result is not None, "Result should not be None" + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + + assert isinstance(circuit_list, list), "circuit_list should be a list" + assert len(circuit_list) == M, f"circuit_list should have {M} elements" + assert isinstance(eu_np_list, list), "eu_np_list should be a list" + assert len(eu_np_list) == M, f"eu_np_list should have {M} elements" + assert isinstance(probs_gptm, np.ndarray), "probs_gptm should be a numpy array" + assert len(probs_gptm) == M, f"probs_gptm should have {M} elements" + assert isinstance(u_choi, np.ndarray), "u_choi should be a numpy array" + assert u_choi.shape == (16, 16), "Choi matrix for 2 qubits should be 16x16" + assert isinstance(u_choi_opt, np.ndarray), "u_choi_opt should be a numpy array" + assert u_choi_opt.shape == ( + 16, + 16, + ), "Optimal Choi matrix for 2 qubits should be 16x16" + assert np.allclose(np.sum(probs_gptm), 1.0), "Probabilities should sum to 1" + assert np.all(probs_gptm >= 0), "All probabilities should be non-negative" + + +def test_mixed_synthesis_sequential_deterministic(): + """Test that mixed_synthesis_sequential returns consistent results with same seed""" + num_qubits = 1 + eps = 1e-3 + M = 16 + unitary = random_su(num_qubits) + + results = [] + for _ in range(3): + result = mixed_synthesis_sequential(unitary, num_qubits, eps, M, seed=123) + assert result is not None + results.append(result) + + # Check that all results are identical + for i in range(1, len(results)): + circuit_list1, eu_np_list1, probs_gptm1, u_choi1, u_choi_opt1 = results[0] + circuit_list2, eu_np_list2, probs_gptm2, u_choi2, u_choi_opt2 = results[i] + + assert len(circuit_list1) == len(circuit_list2) + assert np.allclose( + probs_gptm1, probs_gptm2 + ), "Probabilities should be identical" + assert np.allclose( + u_choi_opt1, u_choi_opt2 + ), "Choi matrices should be identical" + + +def test_mixed_synthesis_parallel_deterministic(): + """Test that mixed_synthesis_parallel returns consistent results with same seed""" + num_qubits = 1 + eps = 1e-3 + M = 16 + unitary = random_su(num_qubits) + + results = [] + for _ in range(3): + result = mixed_synthesis_parallel(unitary, num_qubits, eps, M, seed=123) + assert result is not None + results.append(result) + + # Check that all results are identical + for i in range(1, len(results)): + circuit_list1, eu_np_list1, probs_gptm1, u_choi1, u_choi_opt1 = results[0] + circuit_list2, eu_np_list2, probs_gptm2, u_choi2, u_choi_opt2 = results[i] + + assert len(circuit_list1) == len(circuit_list2) + assert np.allclose( + probs_gptm1, probs_gptm2 + ), "Probabilities should be identical" + assert np.allclose( + u_choi_opt1, u_choi_opt2 + ), "Choi matrices should be identical" + + +def test_mixed_synthesis_sequential_one_qubit(): + """Test mixed_synthesis_sequential with one qubit""" + num_qubits = 1 + eps = 1e-3 + M = 16 + unitary = random_su(num_qubits) + + result = mixed_synthesis_sequential(unitary, num_qubits, eps, M, seed=123) + + assert result is not None + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + + assert len(circuit_list) == M + assert len(eu_np_list) == M + assert len(probs_gptm) == M + assert u_choi_opt.shape == (4, 4), "Choi matrix for 1 qubit should be 4x4" + + +def test_mixed_synthesis_parallel_one_qubit(): + """Test mixed_synthesis_parallel with one qubit""" + num_qubits = 1 + eps = 1e-3 + M = 16 + unitary = random_su(num_qubits) + + result = mixed_synthesis_parallel(unitary, num_qubits, eps, M, seed=123) + + assert result is not None + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + + assert len(circuit_list) == M + assert len(eu_np_list) == M + assert len(probs_gptm) == M + assert u_choi_opt.shape == (4, 4), "Choi matrix for 1 qubit should be 4x4" + + +def test_mixed_synthesis_sequential_with_numpy_array(): + """Test mixed_synthesis_sequential accepts numpy array input""" + num_qubits = 1 + eps = 1e-3 + M = 16 + unitary_mpmath = random_su(num_qubits) + unitary_np = np.array(unitary_mpmath.tolist(), dtype=complex) + + result = mixed_synthesis_sequential(unitary_np, num_qubits, eps, M, seed=123) + + assert result is not None + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + + assert len(circuit_list) == M + assert len(probs_gptm) == M + + +def test_mixed_synthesis_parallel_with_numpy_array(): + """Test mixed_synthesis_parallel accepts numpy array input""" + num_qubits = 1 + eps = 1e-3 + M = 16 + unitary_mpmath = random_su(num_qubits) + unitary_np = np.array(unitary_mpmath.tolist(), dtype=complex) + + result = mixed_synthesis_parallel(unitary_np, num_qubits, eps, M, seed=123) + + assert result is not None + circuit_list, eu_np_list, probs_gptm, u_choi, u_choi_opt = result + + assert len(circuit_list) == M + assert len(probs_gptm) == M + + +if __name__ == "__main__": + pytest.main([__file__, "-v"]) + print("All tests passed successfully!") diff --git a/tests/test_multi_qubit_unitary_approximation.py b/tests/test_multi_qubit_unitary_approximation.py new file mode 100644 index 0000000..8b5075a --- /dev/null +++ b/tests/test_multi_qubit_unitary_approximation.py @@ -0,0 +1,90 @@ +import mpmath +import pytest + +from pygridsynth.config import GridsynthConfig +from pygridsynth.domega_unitary import DOmegaMatrix +from pygridsynth.multi_qubit_unitary_approximation import ( + approximate_multi_qubit_unitary, +) +from pygridsynth.mymath import all_close +from pygridsynth.quantum_circuit import QuantumCircuit + + +def test_approximate_one_qubit_unitary_identity(): + U = mpmath.matrix([[1, 0], [0, 1]]) + num_qubits = 1 + epsilon = "0.0001" + + circuit, approx_unitary = approximate_multi_qubit_unitary(U, num_qubits, epsilon) + + assert isinstance(circuit, QuantumCircuit) + assert isinstance(approx_unitary, mpmath.matrix) + assert len(circuit) == 0 + assert all_close(approx_unitary, U, tol=float(epsilon)) + + +def test_approximate_two_qubit_unitary(): + U = mpmath.matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) + num_qubits = 2 + epsilon = "0.0001" + + circuit, approx_unitary = approximate_multi_qubit_unitary(U, num_qubits, epsilon) + + assert isinstance(circuit, QuantumCircuit) + assert isinstance(approx_unitary, mpmath.matrix) + assert len(circuit) > 0 + assert all_close(approx_unitary, U, tol=float(epsilon)) + + +def test_approximate_three_qubit_unitary_identity(): + num_qubits = 3 + n = 2**num_qubits + U = mpmath.eye(n) + epsilon = "0.0001" + + circuit, approx_unitary = approximate_multi_qubit_unitary( + U, num_qubits, epsilon, return_domega_matrix=True + ) + + assert isinstance(circuit, QuantumCircuit) + assert isinstance(approx_unitary, DOmegaMatrix) + assert all_close(approx_unitary.to_complex_matrix, U, tol=float(epsilon)) + + +def test_approximate_multi_qubit_unitary_return_matrix_type(): + U = mpmath.matrix([[1, 0], [0, 1]]) + num_qubits = 1 + epsilon = "0.0001" + + circuit, result_matrix = approximate_multi_qubit_unitary( + U, num_qubits, epsilon, return_domega_matrix=False + ) + assert isinstance(circuit, QuantumCircuit) + assert isinstance(result_matrix, mpmath.matrix) + + circuit, result_domega = approximate_multi_qubit_unitary( + U, num_qubits, epsilon, return_domega_matrix=True + ) + assert isinstance(circuit, QuantumCircuit) + assert isinstance(result_domega, DOmegaMatrix) + + +def test_approximate_multi_qubit_unitary_config_handling(): + U = mpmath.matrix([[1, 0], [0, 1]]) + num_qubits = 1 + epsilon = "0.0001" + + custom_cfg = GridsynthConfig(dps=50, measure_time=True) + with pytest.warns( + UserWarning, match="When 'cfg' is provided, 'kwargs' are ignored." + ): + approximate_multi_qubit_unitary( + U, num_qubits, epsilon, cfg=custom_cfg, verbose=100 + ) + + approximate_multi_qubit_unitary(U, num_qubits, epsilon, verbose=1) + + +if __name__ == "__main__": + pytest.main([__file__, "-v"]) + print("All tests passed successfully!")