Question about the demo - "Linear combination of unitaries and block encodings"

Demos
1

Hi,

I was testing out the demo titled “Linear combination of unitaries and block encodings” (authored by @Diego @Jay_Soni @Juan Miguel Arrazola) with the following matrix

[[1. 0. -0.25 0.75]

[0. 1. 0.75 0.25]

[-0.25 0.75 2. 0.]

[0.75 0.25 0. 2.]]

note: the matrix is Hermitian.

The LCU decomposition generates negative coefficients, which adversely affect the square root operation in the following step.

  • Is there any way around this? The block encoding stage seems impossible to proceed to with this matrix. Can someone comment on this?

Thanks.

=======
a11 = np.eye(2)
a22 = np.eye(2)*2

a = 0.25
b = 0.75
a12 = np.array([[-a, b],[b,a]])
a21 = np.matrix.getH(a12)

bb2 = np.block([[a11, a12],
[a21, a22 ]])

2

Hi @arund

Thanks for the question and your interest in the demo!

Actually there is a neat property of unitary matrices that they remain unitary under multiplication by complex phases. In practice this means that we can always absorb the negative sign (or in general complex nature) of any coefficient into the associated unitary operator in the LCU decomposition.

I got your example working using the following adjustments:

  • I used PennyLane’s qml.GlobalPhase operator and the qml.prod function to define the new unitary operators with the minus sign absorbed in.
  • This time the LCU had 4 terms instead of the 2 in the demo, this required additional control qubits to prepare the state (4 coefficients = 2 qubits to represent the state).
  • Additionally, one thing we omitted for simplicity in the demo was the fact that when we use this kind of a block-encoding scheme we end up encoding a normalized version of A. Specifically, we encode (A / 1-norm(A)). So when printing the matrix I multiplied by this factor as well.

Here is the code example:

import numpy as np
import pennylane as qml

a = 0.25
b = 0.75

a11 = np.eye(2)
a22 = np.eye(2)*2

a12 = np.array([[-a, b],[b,a]])
a21 = np.matrix.getH(a12)

A = np.block([[a11, a12],
[a21, a22 ]])

A_norm = np.linalg.norm(A, ord=1)
print(A, "\n")
print("1-norm:", A_norm, "\n")
LCU = qml.pauli_decompose(A)
LCU_coeffs, LCU_ops = LCU.terms()

print(f"LCU decomposition:\n {LCU}")
print(f"Coefficients:\n {LCU_coeffs}")
print(f"Unitaries:\n {LCU_ops}")

## Manual Post-processing: 

# Note that a we can always absorb complex phases into our unitary
# i.e. if U is unitary, then U_tilda = e^(i*phi) * U is also unitary. 

new_LCU_coeffs = []
new_LCU_ops = [] 

for coeff, op in zip(LCU_coeffs, LCU_ops):
    
    if coeff < 0: 
        coeff = -1 * coeff
        op = qml.prod(qml.GlobalPhase(np.pi, wires=op.wires), op)  # e^i*pi = -1
    
    new_LCU_coeffs.append(coeff)
    new_LCU_ops.append(op)

print(f"\nNew Coefficients:\n {new_LCU_coeffs}")
print(f"New Unitaries:\n {new_LCU_ops}", "\n \n")

alphas = (np.sqrt(new_LCU_coeffs) / np.linalg.norm(np.sqrt(new_LCU_coeffs)))

# unitaries
ops = new_LCU_ops
# relabeling wires: 0 → 2, and 1 → 3
unitaries = [qml.map_wires(op, {0: 2, 1: 3}) for op in ops]

dev2 = qml.device("default.qubit", wires=[0,1,2,3])
@qml.qnode(dev2)
def lcu_circuit():  # block_encode
    # PREP
    qml.StatePrep(alphas, wires=[0, 1])

    # SEL
    qml.Select(unitaries, control=[0, 1])

    # PREP_dagger
    qml.adjoint(qml.StatePrep(alphas, wires=[0, 1]))
    return qml.state()


output_matrix = qml.matrix(lcu_circuit)()
print("A:\n", A, "\n")
print("Block-encoded A:\n")
print(np.real(np.round((A_norm * output_matrix)[:4, :4],2)))

Let me know if this answers your question,

Thanks,

2 Likes
3

@Jay_Soni Jay, thank you for putting in so much effort. It answers the core of the question. I am, in fact, investigating the technique for a bunch of structured matrices. So the post-processing step is indeed necessary.
In hindsight, I should have posted a longer version of the code. Since I did not want to clutter the question I restrained myself. Indeed, it had more terms in the decomposition. The 1-norm was later clarified from another article. I may have a short follow-up later.
Thanks a lot again for creating the demo. I will post the output below so that anyone following the post can see what to expect.

[[ 1.    0.   -0.25  0.75]
 [ 0.    1.    0.75  0.25]
 [-0.25  0.75  2.    0.  ]
 [ 0.75  0.25  0.    2.  ]] 

1-norm: 3.0 

LCU decomposition:
 1.5 * (I(0) @ I(1)) + 0.75 * (X(0) @ X(1)) + -0.25 * (X(0) @ Z(1)) + -0.5 * (Z(0) @ I(1))
Coefficients:
 [ 1.5   0.75 -0.25 -0.5 ]
Unitaries:
 [I(0) @ I(1), X(0) @ X(1), X(0) @ Z(1), Z(0) @ I(1)]

New Coefficients:
 [np.float64(1.5), np.float64(0.75), np.float64(0.25), np.float64(0.5)]
New Unitaries:
 [I(0) @ I(1), X(0) @ X(1), GlobalPhase(3.141592653589793, wires=[0, 1]) @ (X(0) @ Z(1)), GlobalPhase(3.141592653589793, wires=[0, 1]) @ (Z(0) @ I(1))] 
 

A:
 [[ 1.    0.   -0.25  0.75]
 [ 0.    1.    0.75  0.25]
 [-0.25  0.75  2.    0.  ]
 [ 0.75  0.25  0.    2.  ]] 

Block-encoded A:

[[ 1.    0.   -0.25  0.75]
 [ 0.    1.    0.75  0.25]
 [-0.25  0.75  2.    0.  ]
 [ 0.75  0.25  0.    2.  ]]
2 Likes