Source code for pennylane.ops.op_math.controlled_decompositions
# Copyright 2018-2023 Xanadu Quantum Technologies Inc.# Licensed under the Apache License, Version 2.0 (the "License");# you may not use this file except in compliance with the License.# You may obtain a copy of the License at# http://www.apache.org/licenses/LICENSE-2.0# Unless required by applicable law or agreed to in writing, software# distributed under the License is distributed on an "AS IS" BASIS,# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.# See the License for the specific language governing permissions and# limitations under the License."""This submodule defines functions to decompose controlled operations"""fromcopyimportcopyfromtypingimportOptionalimportnumpyasnpimportnumpy.linalgasnplimportpennylaneasqmlfrompennylaneimportmathfrompennylane.operationimportOperation,Operatorfrompennylane.ops.op_math.decompositions.single_qubit_unitaryimport(_get_single_qubit_rot_angles_via_matrix,)frompennylane.wiresimportWiresdef_is_single_qubit_special_unitary(op):mat=op.matrix()det=mat[0,0]*mat[1,1]-mat[0,1]*mat[1,0]returnqml.math.allclose(det,1)def_convert_to_su2(U,return_global_phase=False):r"""Convert a 2x2 unitary matrix to :math:`SU(2)`. Args: U (array[complex]): A matrix, presumed to be :math:`2 \times 2` and unitary. return_global_phase (bool): If `True`, the return will include the global phase. If `False`, only the :math:`SU(2)` representative is returned. Returns: array[complex]: A :math:`2 \times 2` matrix in :math:`SU(2)` that is equivalent to U up to a global phase. If ``return_global_phase=True``, a 2-element tuple is returned, with the first element being the :math:`SU(2)` equivalent and the second, the global phase. """# Compute the determinantswithnp.errstate(divide="ignore",invalid="ignore"):dets=math.linalg.det(U)global_phase=math.cast_like(math.angle(dets),1.0)/2U_SU2=math.cast_like(U,dets)*math.exp(-1j*global_phase)return(U_SU2,global_phase)ifreturn_global_phaseelseU_SU2def_convert_to_real_diagonal(q:np.ndarray)->np.ndarray:""" Change the phases of Q so the main diagonal is real, and return the modified Q. """exp_angles=np.angle(np.diag(q))returnq*np.exp(-1j*exp_angles).reshape((1,2))def_param_su2(ar:float,ai:float,br:float,bi:float):""" Create a matrix in the SU(2) form from complex parameters a, b. The resulting matrix is not guaranteed to be in SU(2), unless |a|^2 + |b|^2 = 1. """returnnp.array([[ar+1j*ai,-br+1j*bi],[br+1j*bi,ar+1j*-ai]])def_bisect_compute_a(u:np.ndarray):""" Given the U matrix, compute the A matrix such that At x A x At x A x = U where At is the adjoint of A and x is the Pauli X matrix. """x=np.real(u[0,1])z=u[1,1]zr=np.real(z)zi=np.imag(z)ifnp.isclose(zr,-1):# special case [[-1, 0], [0, -1]]# would cause divide by 0 with the other formula, so we use hardcoded solutionreturnnp.array([[1,-1],[1,1]])*2**-0.5ar=np.sqrt((np.sqrt((zr+1)/2)+1)/2)mul=1/(2*np.sqrt((zr+1)*(np.sqrt((zr+1)/2)+1)))ai=zi*mulbr=x*mulbi=0return_param_su2(ar,ai,br,bi)def_bisect_compute_b(u:np.ndarray):""" Given the U matrix, compute the B matrix such that H Bt x B x H = U where Bt is the adjoint of B, H is the Hadamard matrix, and x is the Pauli X matrix. """sqrt=np.sqrtAbs=np.absw=np.real(u[0,0])s=np.real(u[1,0])t=np.imag(u[1,0])ifnp.isclose(s,0):b=0ifnp.isclose(t,0):ifw<0:c=0d=sqrt(-w)else:c=sqrt(w)d=0else:c=sqrt(2-2*w)*(-w/2-1/2)/td=sqrt(2-2*w)/2elifnp.isclose(t,0):b=(1/2-w/2)*sqrt(2*w+2)/sc=sqrt(2*w+2)/2d=0else:b=sqrt(2)*s*sqrt((1-w)/(s**2+t**2))*Abs(t)/(2*t)c=sqrt(2)*sqrt((1-w)/(s**2+t**2))*(w+1)*Abs(t)/(2*t)d=-sqrt(2)*sqrt((1-w)/(s**2+t**2))*Abs(t)/2return_param_su2(c,d,b,0)def_multi_controlled_zyz(rot_angles,global_phase,target_wire:Wires,control_wires:Wires,work_wires:Optional[Wires]=None,)->list[Operator]:# The decomposition of zyz for special unitaries with multiple control wires# defined in Lemma 7.9 of https://arxiv.org/pdf/quant-ph/9503016ifnotqml.math.allclose(0.0,global_phase,atol=1e-6,rtol=0):raiseValueError(f"The global_phase should be zero, instead got: {global_phase}.")# Unpack the rotation anglesphi,theta,omega=rot_angles# We use the conditional statements to account when decomposition is ran within a queuedecomp=[]cop_wires=(control_wires[-1],target_wire[0])# Add operator Aifnotqml.math.allclose(0.0,phi,atol=1e-8,rtol=0):decomp.append(qml.CRZ(phi,wires=cop_wires))ifnotqml.math.allclose(0.0,theta/2,atol=1e-8,rtol=0):decomp.append(qml.CRY(theta/2,wires=cop_wires))decomp.append(qml.ctrl(qml.X(target_wire),control=control_wires[:-1],work_wires=work_wires))# Add operator Bifnotqml.math.allclose(0.0,theta/2,atol=1e-8,rtol=0):decomp.append(qml.CRY(-theta/2,wires=cop_wires))ifnotqml.math.allclose(0.0,-(phi+omega)/2,atol=1e-6,rtol=0):decomp.append(qml.CRZ(-(phi+omega)/2,wires=cop_wires))decomp.append(qml.ctrl(qml.X(target_wire),control=control_wires[:-1],work_wires=work_wires))# Add operator Cifnotqml.math.allclose(0.0,(omega-phi)/2,atol=1e-8,rtol=0):decomp.append(qml.CRZ((omega-phi)/2,wires=cop_wires))returndecompdef_single_control_zyz(rot_angles,global_phase,target_wire,control_wires:Wires):# The zyz decomposition of a general unitary with single control wire# defined in Lemma 7.9 of https://arxiv.org/pdf/quant-ph/9503016# Unpack the rotation anglesphi,theta,omega=rot_angles# We use the conditional statements to account when decomposition is ran within a queuedecomp=[]# Add negative of global phase. Compare definition of qml.GlobalPhase and Ph(delta) from section 4.1 of Barenco et al.ifnotqml.math.allclose(0.0,global_phase,atol=1e-8,rtol=0):decomp.append(qml.ctrl(qml.GlobalPhase(phi=-global_phase,wires=target_wire),control=control_wires))# Add operator Aifnotqml.math.allclose(0.0,phi,atol=1e-8,rtol=0):decomp.append(qml.RZ(phi,wires=target_wire))ifnotqml.math.allclose(0.0,theta/2,atol=1e-8,rtol=0):decomp.append(qml.RY(theta/2,wires=target_wire))decomp.append(qml.ctrl(qml.X(target_wire),control=control_wires))# Add operator Bifnotqml.math.allclose(0.0,theta/2,atol=1e-8,rtol=0):decomp.append(qml.RY(-theta/2,wires=target_wire))ifnotqml.math.allclose(0.0,-(phi+omega)/2,atol=1e-6,rtol=0):decomp.append(qml.RZ(-(phi+omega)/2,wires=target_wire))decomp.append(qml.ctrl(qml.X(target_wire),control=control_wires))# Add operator Cifnotqml.math.allclose(0.0,(omega-phi)/2,atol=1e-8,rtol=0):decomp.append(qml.RZ((omega-phi)/2,wires=target_wire))returndecomp
[docs]defctrl_decomp_zyz(target_operation:Operator,control_wires:Wires,work_wires:Optional[Wires]=None)->list[Operator]:"""Decompose the controlled version of a target single-qubit operation This function decomposes both single and multiple controlled single-qubit target operations using the decomposition defined in Lemma 4.3 and Lemma 5.1 for single `controlled_wires`, and Lemma 7.9 for multiple `controlled_wires` from `Barenco et al. (1995) <https://arxiv.org/abs/quant-ph/9503016>`_. Args: target_operation (~.operation.Operator): the target operation to decompose control_wires (~.wires.Wires): the control wires of the operation. Returns: list[Operation]: the decomposed operations Raises: ValueError: if ``target_operation`` is not a single-qubit operation **Example** We can create a controlled operation using `qml.ctrl`, or by creating the decomposed controlled version of using `qml.ctrl_decomp_zyz`. .. code-block:: python import pennylane as qml dev = qml.device("default.qubit", wires=2) @qml.qnode(dev) def expected_circuit(op): qml.Hadamard(wires=0) qml.ctrl(op, [0]) return qml.probs() @qml.qnode(dev) def decomp_circuit(op): qml.Hadamard(wires=0) qml.ops.ctrl_decomp_zyz(op, [0]) return qml.probs() Measurements on both circuits will give us the same results: >>> op = qml.RX(0.123, wires=1) >>> expected_circuit(op) tensor([0.5 , 0. , 0.49811126, 0.00188874], requires_grad=True) >>> decomp_circuit(op) tensor([0.5 , 0. , 0.49811126, 0.00188874], requires_grad=True) """iflen(target_operation.wires)!=1:raiseValueError("The target operation must be a single-qubit operation, instead "f"got {target_operation.__class__.__name__}.")control_wires=Wires(control_wires)target_wire=target_operation.wiresifisinstance(target_operation,Operation):try:rot_angles=target_operation.single_qubit_rot_angles()exceptNotImplementedError:rot_angles=_get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))else:rot_angles=_get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))_,global_phase=_convert_to_su2(qml.matrix(target_operation),return_global_phase=True)return(_multi_controlled_zyz(rot_angles,global_phase,target_wire,control_wires,work_wires)iflen(control_wires)>1else_single_control_zyz(rot_angles,global_phase,target_wire,control_wires))
def_ctrl_decomp_bisect_od(u:np.ndarray,target_wire:Wires,control_wires:Wires,):"""Decompose the controlled version of a target single-qubit operation Not backpropagation compatible (as currently implemented). Use only with numpy. This function decomposes a controlled single-qubit target operation using the decomposition defined in section 3.1, Theorem 1 of `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_. The target operation's matrix must have a real off-diagonal for this specialized method to work. .. warning:: This method will add a global phase for target operations that do not belong to the SU(2) group. Args: u (np.ndarray): the target operation's matrix target_wire (~.wires.Wires): the target wire of the controlled operation control_wires (~.wires.Wires): the control wires of the operation Returns: list[Operation]: the decomposed operations Raises: ValueError: if ``u`` does not have a real off-diagonal """ui=np.imag(u)ifnotnp.isclose(ui[1,0],0)ornotnp.isclose(ui[0,1],0):raiseValueError(f"Target operation's matrix must have real off-diagonal, but it is {u}")a=_bisect_compute_a(u)mid=(len(control_wires)+1)//2# for odd n, make control_k1 biggercontrol_k1=control_wires[:mid]control_k2=control_wires[mid:]defcomponent():return[qml.ctrl(qml.X(target_wire),control=control_k1,work_wires=control_k2),qml.QubitUnitary(a,target_wire),qml.ctrl(qml.X(target_wire),control=control_k2,work_wires=control_k1),qml.adjoint(qml.QubitUnitary(a,target_wire)),]returncomponent()+component()def_ctrl_decomp_bisect_md(u:np.ndarray,target_wire:Wires,control_wires:Wires,):"""Decompose the controlled version of a target single-qubit operation Not backpropagation compatible (as currently implemented). Use only with numpy. This function decomposes a controlled single-qubit target operation using the decomposition defined in section 3.1, Theorem 2 of `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_. The target operation's matrix must have a real main-diagonal for this specialized method to work. .. warning:: This method will add a global phase for target operations that do not belong to the SU(2) group. Args: u (np.ndarray): the target operation's matrix target_wire (~.wires.Wires): the target wire of the controlled operation control_wires (~.wires.Wires): the control wires of the operation Returns: list[Operation]: the decomposed operations Raises: ValueError: if ``u`` does not have a real main-diagonal """ui=np.imag(u)ifnotnp.isclose(ui[0,0],0)ornotnp.isclose(ui[1,1],0):raiseValueError(f"Target operation's matrix must have real main-diagonal, but it is {u}")h_matrix=qml.Hadamard.compute_matrix()mod_u=h_matrix@u@h_matrixdecomposition=[qml.Hadamard(target_wire)]decomposition+=_ctrl_decomp_bisect_od(mod_u,target_wire,control_wires)decomposition.append(qml.Hadamard(target_wire))returndecompositiondef_ctrl_decomp_bisect_general(u:np.ndarray,target_wire:Wires,control_wires:Wires,):"""Decompose the controlled version of a target single-qubit operation Not backpropagation compatible (as currently implemented). Use only with numpy. This function decomposes a controlled single-qubit target operation using the decomposition defined in section 3.2 of `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_. .. warning:: This method will add a global phase for target operations that do not belong to the SU(2) group. Args: u (np.ndarray): the target operation's matrix target_wire (~.wires.Wires): the target wire of the controlled operation control_wires (~.wires.Wires): the control wires of the operation Returns: list[Operation]: the decomposed operations """x_matrix=qml.X.compute_matrix()h_matrix=qml.Hadamard.compute_matrix()alternate_h_matrix=x_matrix@h_matrix@x_matrixd,q=npl.eig(u)d=np.diag(d)q=_convert_to_real_diagonal(q)b=_bisect_compute_b(q)c1=b@alternate_h_matrixc2t=b@h_matrixmid=(len(control_wires)+1)//2# for odd n, make control_k1 biggercontrol_k1=control_wires[:mid]control_k2=control_wires[mid:]component=[qml.QubitUnitary(c2t,target_wire),qml.ctrl(qml.X(target_wire),control=control_k2,work_wires=control_k1),qml.adjoint(qml.QubitUnitary(c1,target_wire)),qml.ctrl(qml.X(target_wire),control=control_k1,work_wires=control_k2),]od_decomp=_ctrl_decomp_bisect_od(d,target_wire,control_wires)# cancel two identical multicontrolled x gatesqml.QueuingManager.remove(component[3])qml.QueuingManager.remove(od_decomp[0])adjoint_component=[qml.adjoint(copy(op),lazy=False)foropinreversed(component)]returncomponent[0:3]+od_decomp[1:]+adjoint_component
[docs]defctrl_decomp_bisect(target_operation:Operator,control_wires:Wires,):"""Decompose the controlled version of a target single-qubit operation Not backpropagation compatible (as currently implemented). Use only with numpy. Automatically selects the best algorithm based on the matrix (uses specialized more efficient algorithms if the matrix has a certain form, otherwise falls back to the general algorithm). These algorithms are defined in section 3.1 and 3.2 of `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_. .. warning:: This method will add a global phase for target operations that do not belong to the SU(2) group. Args: target_operation (~.operation.Operator): the target operation to decompose control_wires (~.wires.Wires): the control wires of the operation Returns: list[Operation]: the decomposed operations Raises: ValueError: if ``target_operation`` is not a single-qubit operation **Example:** >>> op = qml.T(0) # uses OD algorithm >>> print(qml.draw(ctrl_decomp_bisect, wire_order=(0,1,2,3,4,5), show_matrices=False)(op, (1,2,3,4,5))) 0: ─╭X──U(M0)─╭X──U(M0)†─╭X──U(M0)─╭X──U(M0)†─┤ 1: ─├●────────│──────────├●────────│──────────┤ 2: ─├●────────│──────────├●────────│──────────┤ 3: ─╰●────────│──────────╰●────────│──────────┤ 4: ───────────├●───────────────────├●─────────┤ 5: ───────────╰●───────────────────╰●─────────┤ >>> op = qml.QubitUnitary([[0,1j],[1j,0]], 0) # uses MD algorithm >>> print(qml.draw(ctrl_decomp_bisect, wire_order=(0,1,2,3,4,5), show_matrices=False)(op, (1,2,3,4,5))) 0: ──H─╭X──U(M0)─╭X──U(M0)†─╭X──U(M0)─╭X──U(M0)†──H─┤ 1: ────├●────────│──────────├●────────│─────────────┤ 2: ────├●────────│──────────├●────────│─────────────┤ 3: ────╰●────────│──────────╰●────────│─────────────┤ 4: ──────────────├●───────────────────├●────────────┤ 5: ──────────────╰●───────────────────╰●────────────┤ >>> op = qml.Hadamard(0) # uses general algorithm >>> print(qml.draw(ctrl_decomp_bisect, wire_order=(0,1,2,3,4,5), show_matrices=False)(op, (1,2,3,4,5))) 0: ──U(M0)─╭X──U(M1)†──U(M2)─╭X──U(M2)†─╭X──U(M2)─╭X──U(M2)†─╭X──U(M1)─╭X──U(M0)─┤ 1: ────────│─────────────────│──────────├●────────│──────────├●────────│─────────┤ 2: ────────│─────────────────│──────────├●────────│──────────├●────────│─────────┤ 3: ────────│─────────────────│──────────╰●────────│──────────╰●────────│─────────┤ 4: ────────├●────────────────├●───────────────────├●───────────────────├●────────┤ 5: ────────╰●────────────────╰●───────────────────╰●───────────────────╰●────────┤ """iflen(target_operation.wires)>1:raiseValueError("The target operation must be a single-qubit operation, instead "f"got {target_operation}.")target_matrix=target_operation.matrix()target_wire=target_operation.wirestarget_matrix=_convert_to_su2(target_matrix)target_matrix_imag=np.imag(target_matrix)ifnp.isclose(target_matrix_imag[1,0],0)andnp.isclose(target_matrix_imag[0,1],0):# Real off-diagonal specialized algorithm - 16n+O(1) CNOTsreturn_ctrl_decomp_bisect_od(target_matrix,target_wire,control_wires)ifnp.isclose(target_matrix_imag[0,0],0)andnp.isclose(target_matrix_imag[1,1],0):# Real main-diagonal specialized algorithm - 16n+O(1) CNOTsreturn_ctrl_decomp_bisect_md(target_matrix,target_wire,control_wires)# General algorithm - 20n+O(1) CNOTsreturn_ctrl_decomp_bisect_general(target_matrix,target_wire,control_wires)
defdecompose_mcx(control_wires,target_wire,work_wires):"""Decomposes the multi-controlled PauliX gate using decompositions from `Barenco et al. (1995) <https://arxiv.org/abs/quant-ph/9503016>`_"""num_work_wires_needed=len(control_wires)-2iflen(control_wires)==1:return[qml.CNOT(wires=control_wires+Wires(target_wire))]iflen(control_wires)==2:returnqml.Toffoli.compute_decomposition(wires=control_wires+Wires(target_wire))iflen(work_wires)>=num_work_wires_needed:# Lemma 7.2return_decompose_mcx_with_many_workers(control_wires,target_wire,work_wires)iflen(work_wires)>=1:# Lemma 7.3return_decompose_mcx_with_one_worker(control_wires,target_wire,work_wires[0])# Lemma 7.5withqml.QueuingManager.stop_recording():op=qml.X(target_wire)return_decompose_multicontrolled_unitary(op,control_wires)def_decompose_multicontrolled_unitary(op,control_wires):"""Decomposes general multi controlled unitary with no work wires Follows approach from Lemma 7.5 combined with 7.3 and 7.2 of https://arxiv.org/abs/quant-ph/9503016. We are assuming this decomposition is used only in the general cases """ifnotop.has_matrixorlen(op.wires)!=1:raiseValueError("The target operation must be a single-qubit operation with a matrix representation")target_wire=op.wiresiflen(control_wires)==0:return[op]iflen(control_wires)==1:returnctrl_decomp_zyz(op,control_wires)if_is_single_qubit_special_unitary(op):returnctrl_decomp_bisect(op,control_wires)# use recursive decomposition of general gatereturn_decompose_recursive(op,1.0,control_wires,target_wire,Wires([]))def_decompose_recursive(op,power,control_wires,target_wire,work_wires):"""Decompose multicontrolled operator recursively using lemma 7.5 Number of gates in decomposition are: O(len(control_wires)^2) """iflen(control_wires)==1:withqml.QueuingManager.stop_recording():powered_op=qml.pow(op,power,lazy=True)returnctrl_decomp_zyz(powered_op,control_wires)withqml.QueuingManager.stop_recording():cnots=decompose_mcx(control_wires=control_wires[:-1],target_wire=control_wires[-1],work_wires=work_wires+target_wire,)withqml.QueuingManager.stop_recording():powered_op=qml.pow(op,0.5*power,lazy=True)powered_op_adj=qml.adjoint(powered_op,lazy=True)ifqml.QueuingManager.recording():decomposition=[*ctrl_decomp_zyz(powered_op,control_wires[-1]),*(qml.apply(o)foroincnots),*ctrl_decomp_zyz(powered_op_adj,control_wires[-1]),*(qml.apply(o)foroincnots),*_decompose_recursive(op,0.5*power,control_wires[:-1],target_wire,control_wires[-1]+work_wires),]else:decomposition=[*ctrl_decomp_zyz(powered_op,control_wires[-1]),*cnots,*ctrl_decomp_zyz(powered_op_adj,control_wires[-1]),*cnots,*_decompose_recursive(op,0.5*power,control_wires[:-1],target_wire,control_wires[-1]+work_wires),]returndecompositiondef_decompose_mcx_with_many_workers(control_wires,target_wire,work_wires):"""Decomposes the multi-controlled PauliX gate using the approach in Lemma 7.2 of https://arxiv.org/abs/quant-ph/9503016, which requires a suitably large register of work wires"""num_work_wires_needed=len(control_wires)-2work_wires=work_wires[:num_work_wires_needed]work_wires_reversed=list(reversed(work_wires))control_wires_reversed=list(reversed(control_wires))gates=[]foriinrange(len(work_wires)):ctrl1=control_wires_reversed[i]ctrl2=work_wires_reversed[i]t=target_wireifi==0elsework_wires_reversed[i-1]gates.append(qml.Toffoli(wires=[ctrl1,ctrl2,t]))gates.append(qml.Toffoli(wires=[*control_wires[:2],work_wires[0]]))foriinreversed(range(len(work_wires))):ctrl1=control_wires_reversed[i]ctrl2=work_wires_reversed[i]t=target_wireifi==0elsework_wires_reversed[i-1]gates.append(qml.Toffoli(wires=[ctrl1,ctrl2,t]))foriinrange(len(work_wires)-1):ctrl1=control_wires_reversed[i+1]ctrl2=work_wires_reversed[i+1]t=work_wires_reversed[i]gates.append(qml.Toffoli(wires=[ctrl1,ctrl2,t]))gates.append(qml.Toffoli(wires=[*control_wires[:2],work_wires[0]]))foriinreversed(range(len(work_wires)-1)):ctrl1=control_wires_reversed[i+1]ctrl2=work_wires_reversed[i+1]t=work_wires_reversed[i]gates.append(qml.Toffoli(wires=[ctrl1,ctrl2,t]))returngatesdef_decompose_mcx_with_one_worker(control_wires,target_wire,work_wire):"""Decomposes the multi-controlled PauliX gate using the approach in Lemma 7.3 of https://arxiv.org/abs/quant-ph/9503016, which requires a single work wire"""tot_wires=len(control_wires)+2partition=int(np.ceil(tot_wires/2))first_part=control_wires[:partition]second_part=control_wires[partition:]gates=[qml.ctrl(qml.X(work_wire),control=first_part,work_wires=second_part+target_wire),qml.ctrl(qml.X(target_wire),control=second_part+work_wire,work_wires=first_part),qml.ctrl(qml.X(work_wire),control=first_part,work_wires=second_part+target_wire),qml.ctrl(qml.X(target_wire),control=second_part+work_wire,work_wires=first_part),]returngates