Consider the quantum operation of addition mod M:
(x,y)->(x,(x+y)%M), where x and y
are m-qubit pieces
of an n-qubit register, with n=2m,
N=2n, and M=2m.
A high-level simulation can be performed simply using classical addition. In C pseudo-code:
I = state index
x = I >> m; // top m bits
y = I & mask; // bottom m bits, mask = 2m-1 = 0111..1 (m 1's)
j = (x + y) % M; // classical addition
J = (x << m) | j; // concatenate the bits
tj = a[J]; // temporarily save original a[J]
a[J] = a[I]; // update a[J]
where a is the array of quantum amplitudes.
The amplitude of state I=(x,y) is copied
to state J=(x,(x+y)%M), while storing the original
state J amplitude in a temporary variable tj.
This code will be in a loop going through the state indices.
To avoid using an array of temporary variables, tj
must be copied to the appropriate location in the next iteration
before being used again to save the next a[J] amplitude.
The following shows how this is done in the actual C++ code
add.cc.
Running the program without any options will initialize the state amplitudes to the
(unnormalized) value of the sum. For example, with m=1,
the amplitude for state (1,1) is initialized to (1,0)=2
since it will be moved to state (1,(1+1)%2)=(1,0) after performing the sum.
$ ./add 1
before, a[i] = sum: x y x j 0 + 0 -> 0 0 : 0 0 + 1 -> 0 1 : 1 1 + 0 -> 1 1 : 3 1 + 1 -> 1 0 : 2 |
after, a[i] = i: x j 0 0 : 0 0 1 : 1 1 0 : 2 1 1 : 3 |
a[i] should
equal i and the amplitudes should be in sequence 0, 1, ..., N-1.
$ ./add 2
before, a[i] = sum: x y x j 00 + 00 -> 00 00 : 0 00 + 01 -> 00 01 : 1 00 + 10 -> 00 10 : 2 00 + 11 -> 00 11 : 3 01 + 00 -> 01 01 : 5 01 + 01 -> 01 10 : 6 01 + 10 -> 01 11 : 7 01 + 11 -> 01 00 : 4 10 + 00 -> 10 10 : 10 10 + 01 -> 10 11 : 11 10 + 10 -> 10 00 : 8 10 + 11 -> 10 01 : 9 11 + 00 -> 11 11 : 15 11 + 01 -> 11 00 : 12 11 + 10 -> 11 01 : 13 11 + 11 -> 11 10 : 14 |
after, a[i] = i: x j 00 00 : 0 00 01 : 1 00 10 : 2 00 11 : 3 01 00 : 4 01 01 : 5 01 10 : 6 01 11 : 7 10 00 : 8 10 01 : 9 10 10 : 10 10 11 : 11 11 00 : 12 11 01 : 13 11 10 : 14 11 11 : 15 |
The before outputs show that for each x,
(x+y)%M is a circular shift of y by amount x.
For example, with M=8,
look at the binary values of y and j in the
./add 3 output below.
The amplitudes are moved
to position j
as follows for x = 1, 2, 4:
x = 1, j = 1, 2, 3, 4, 5, 6, 7, 0 x = 2, j = 2, 4, 6, 0, 3, 5, 7, 1 x = 4, j = 4, 0, 5, 1, 6, 2, 7, 3So for
x=1:
y=0 moves to j=1,
y=1 moves to j=2, etc.
For x=2:
y=0 moves to j=2,
y=2 moves to j=4,
y=4 moves to j=6,
y=6 moves to j=0,
y=1 moves to j=3, etc.
Note the transition: after moving y=6 the process starts over with y=1.
This transition happens more often for x=4:
y=0 moves to j=4,
y=4 moves to j=0,
y=1 moves to j=5,
y=5 moves to j=1,
y=2 moves to j=6, etc.
$ ./add 3
before, a[i] = sum: ... x y x j 001 + 000 -> 001 001 : 9 001 + 001 -> 001 010 : 10 001 + 010 -> 001 011 : 11 001 + 011 -> 001 100 : 12 001 + 100 -> 001 101 : 13 001 + 101 -> 001 110 : 14 001 + 110 -> 001 111 : 15 001 + 111 -> 001 000 : 8 010 + 000 -> 010 010 : 18 010 + 001 -> 010 011 : 19 010 + 010 -> 010 100 : 20 010 + 011 -> 010 101 : 21 010 + 100 -> 010 110 : 22 010 + 101 -> 010 111 : 23 010 + 110 -> 010 000 : 16 010 + 111 -> 010 001 : 17 ... 100 + 000 -> 100 100 : 36 100 + 001 -> 100 101 : 37 100 + 010 -> 100 110 : 38 100 + 011 -> 100 111 : 39 100 + 100 -> 100 000 : 32 100 + 101 -> 100 001 : 33 100 + 110 -> 100 010 : 34 100 + 111 -> 100 011 : 35 ... |
after, a[i] = i: ... x j 001 000 : 8 001 001 : 9 001 010 : 10 001 011 : 11 001 100 : 12 001 101 : 13 001 110 : 14 001 111 : 15 010 000 : 16 010 001 : 17 010 010 : 18 010 011 : 19 010 100 : 20 010 101 : 21 010 110 : 22 010 111 : 23 ... 100 000 : 32 100 001 : 33 100 010 : 34 100 011 : 35 100 100 : 36 100 101 : 37 100 110 : 38 100 111 : 39 ... |
The addition loop starts with x=1 since addition of x=0
has no effect. A two-point sliding window of temporary values ti,tj
is used:
unsigned long M = 1LU << m, mask = M-1; // mask = 0111..1 (m 1's)
for( unsigned long x = 1; x < M; ++x) { // skip x=0
unsigned long start = 0, i = start, j = (i+x) & mask, count = M, xm = x << m;
Complex ti = q.a[xm], tj;
Index i starts at start, which is 0 initially.
But if index j wraps around back to start,
then start will be incremented to handle the transition noted above.
while( 1) { // until count=0
unsigned long xj = xm | j;
tj = q.a[xj]; q.a[xj] = ti; ti = tj; --count; if( count == 0) break;
The current amplitude q.a[xj] is copied to tj and
then overwritten with the temporary value ti from the previous iteration.
ti=tj then slides the window.
if( j != start) i = j; else { ++start; i = start; ti = q.a[xm|i]; }
j = (i+x) & mask;
}
}
Normally i=j moves the index window for the next iteration. But if
index j has wrapped around, then start is incremented
and i and ti are reset.
If a quantum algorithm does not involve entanglement
then the state vector is separable and can be represented using 2n
amplitudes instead of 2n. The simulation of such a
product state would be efficient and practical for large n and there
would be no need to use an actual quantum computer for the computation.
Considering addition, if that was the only computation needed, then it could be done classically or in a simulation just as efficiently as on a quantum computer, if there is no entanglement. Of course addition is generally just one step in an overall algorithm that involves many other types of computations, such as multiplication, modular exponentiation, Fourier transform, etc.
But this raises the question of whether addition itself may create entanglement,
if there was no entanglement present before the addition. To examine this issue,
the add.cc simulator -E option can be use to see the trace distance
entanglement measure for each qubit, where 0 indicates no entanglement and 1 is the maximum
entanglement(1).
In the first example shown above, with state amplitudes initialized to the sum values,
there is entanglement before performing the addition, with E=0.2
for qubit #0 and E=0.692308 for qubit #1:
$ ./add -E 1 before, a[i] = sum: 0 + 0 -> 0 0 : 0 0 + 1 -> 0 1 : 1 1 + 0 -> 1 1 : 3 1 + 1 -> 1 0 : 2 E: 0.2 0.692308 after, a[i] = i: 0 0 : 0 0 1 : 1 1 0 : 2 1 1 : 3 E: 0.1 0.307692After performing the addition, the entanglement measure is
E=0.1 for qubit #0 and E=0.307692 for qubit #1.
If the state amplitudes are initialized to a uniform superposition then there is no entanglement before the addition, and also none after, since permuting the equal amplitudes has no effect:
$ ./add -uE 1
before:
00 0.25 (0.5)
01 0.25 (0.5)
10 0.25 (0.5)
11 0.25 (0.5)
E: 0 0
after:
0 0 : 0.5
0 1 : 0.5
1 0 : 0.5
1 1 : 0.5
E: 0 0
However, if one CNOT is performed after initialization, then there is
maximum entanglement before the addition, and also after:
$ ./add -uCE 1
before:
00 0.25 (0.5)
01 0.25 (0.5)
10 0.25 (0.5)
11 0.25 (-0.5)
E: 1 1
after:
0 0 : 0.5
0 1 : 0.5
1 0 : -0.5
1 1 : 0.5
E: 1 1
To demonstrate that addition can create entanglement, start with a product state that has mostly distinct unnormalized amplitudes. There is no entanglement before the addition, but after the addition all of the qubits are entangled:
$ ./add -pE 1
before:
00 9 (3)
01 36 (6)
10 16 (4)
11 64 (8)
E: 0 0
after:
0 0 : 3
0 1 : 6
1 0 : 8
1 1 : 4
E: 0.341412 0.36
Some more examples, increasing the number of qubits:
$ ./add -pE 1 | egrep '^(before|E:|after)' before: E: 0 0 after: E: 0.341412 0.36 $ ./add -pE 2 | egrep '^(before|E:|after)' before: E: 0 0 0 0 after: E: 0.382632 0.078396 0.385344 0.0784 $ ./add -pE 3 | egrep '^(before|E:|after)' before: E: 0 0 0 0 0 5.96046e-08 after: E: 0.390365 0.093587 0.0325131 0.390269 0.0934461 0.0325179 $ ./add -pE 4 | egrep '^(before|E:|after)' before: E: 0 5.96046e-07 0 0 0 0 2.98023e-07 1.78814e-07 after: E: 0.392432 0.0976034 0.0409784 0.0176182 0.391574 0.0974197 0.0410607 0.0176209 $ ./add -pE 5 | egrep '^(before|E:|after)' before: E: -nan -nan -nan -nan -nan -nan -nan -nan -nan -nan after: E: -nan -nan -nan -nan -nan -nan -nan -nan -nan -nanWith
m=5 and float precision,
the numerator and denominator of the entanglement
computation overflow. The overflow can be avoided by normalizing the amplitudes:
$ ./add -pnE 5 | egrep '^(before|E:|after)' before: E: 0 5.36442e-07 0 1.78814e-07 0 5.96046e-08 0 1.01328e-06 4.17233e-07 7.15256e-07 after: E: 0.393069 0.0988368 0.0436134 0.0229837 0.0110177 0.391977 0.0986472 0.0436966 0.0230408 0.0110207
In [1] a quantum algorithm for SVD is presented, along with applications for state swapping and encoding. Here we will use the implementation from [2] to directly compute the SVD, and demonstrate those applications.
The SVD is performed on a matrix A
representing a partition of the state vector (x,y)
into subsystems x and y. The rows of the matrix correspond to
subsystem x and the columns to subsystem y.
For quantum addition the subsystems are naturally the same size but in general the
partitions could be different sizes.
The SVD decomposition used here produces A=Q'*S*V where
Q is unitary, the rows of V are orthonormal,
and S is diagonal {si} containing the real
non-negative singular values in decreasing order.
The number of non-zero singular values is the rank.
If the rank is 1 then the system is in a product state with no entanglement,
otherwise there is entanglement.
The matrix A can also be written in summation form as
A=∑i=1..rank{si*qi'*vi}
where qi and vi
are rows of Q and V.
The sets {qi} and {vi} are
orthonormal bases for the subsystems x and y.
This is the Schmidt form corresponding to equation (2) in [1].
Consider 2-bit addition with an initial unnormalized product state, and compare the system state representation in vector and matrix forms:
$ ./add -pE 2 before: x y 0000 1.102e+04 (105) 0001 4.41e+04 (210) 0010 1.96e+04 (140) 0011 7.84e+04 (280) 0100 1.588e+04 (126) 0101 6.35e+04 (252) 0110 2.822e+04 (168) 0111 1.129e+05 (336) 1000 1.44e+04 (120) 1001 5.76e+04 (240) 1010 2.56e+04 (160) 1011 1.024e+05 (320) 1100 2.074e+04 (144) 1101 8.294e+04 (288) 1110 3.686e+04 (192) 1111 1.475e+05 (384) E: 0 0 0 0 after: x y 00 00 : 105 00 01 : 210 00 10 : 140 00 11 : 280 01 00 : 336 01 01 : 126 01 10 : 252 01 11 : 168 10 00 : 160 10 01 : 320 10 10 : 120 10 11 : 240 11 00 : 288 11 01 : 192 11 10 : 384 11 11 : 144 E: 0.382632 0.078396 0.385344 0.0784 |
$ ./add -pSS 2 (annotated, Q not shown) # before: y 00 01 10 11 x A = [ 00 105 210 140 280 01 126 252 168 336 10 120 240 160 320 11 144 288 192 384 ]; V = [ 0.268328 0.536656 0.357771 0.715542 -4.88281e-06 -9.76562e-06 3.66211e-06 7.32422e-06 0 0 0 0 -7.32422e-06 -1.46484e-05 5.49316e-06 1.09863e-05 ]; S = [ 928.238 1.86265e-10 0 4.19095e-10 ]; # after: A = [ 105 210 140 280 336 126 252 168 160 320 120 240 288 192 384 144 ]; V = [ 0.527262 0.474151 0.537544 0.456314 -0.433077 0.542525 -0.496307 0.521336 -0.65303 0.329958 0.609177 -0.305913 0.328622 0.609897 -0.305988 -0.652996 ]; S = [ 873.606 285.529 99.8982 83.2497 ]; |
After performing the addition the system is fully entangled and the rank is 4.
The swap described in [1] uses the SVD to produce the transpose
of the state A by operating on each subsystem with the matrix
formed from the basis of the other subsystem. Since the subsystems
x and y are associated with the rows and columns of
A respectively, the transpose represents a state with the
subsystems swapped.
In matrix form, A=Q'*S*V and what [1] refers to
as "applying the QSVD" is equivalent to computing D=Q*A*V',
where D is diagonal and equals S.
Then the swap is performed by computing V'*D'*Q,
which will equal A'.
Example using double precision:
$ ./add -pnSSa 2 > jj; octave jj
>> A
A =
0.11312 0.22624 0.15082 0.30165
0.36198 0.13574 0.27148 0.18099
0.17237 0.34474 0.12928 0.25855
0.31027 0.20684 0.41369 0.15513
>> S
S =
0.94115 0.3076 0.10762 0.089686
>> D = Q*A*V'; D(abs(D)<10*eps)=0
D =
0.94115 0 0 0
0 0.3076 0 0
0 0 0.10762 0
0 0 0 0.089686
>> V'*D'*Q # equals A'
ans =
0.11312 0.36198 0.17237 0.31027
0.22624 0.13574 0.34474 0.20684
0.15082 0.27148 0.12928 0.41369
0.30165 0.18099 0.25855 0.15513
After applying the QSVD, the state matrix D is diagonal which
means that the state of subsystem x is identical to the state of
subsystem y. In the example above
for instance the state amplitudes using binary index notation are:
a[00,00]=0.94115, a[01,01]=0.3076, a[10,10]=0.10762, a[11,11]=0.089686.
This means that one could "throw away" one of the subsystems without losing any
information. However, since it is not possible to clone an arbitrary quantum
state [3], the original bipartite quantum state can
not be recovered from the remaining subsystem.
The no-cloning restriction is overcome using the quantum encoder described in
[1]. By applying a CNOT operation between each pair of
coincident qubits in the subsystems, controlled at subsystem x
and targeted to subsystem y for example, the state of subsystem y
becomes 0. So then subsystem y could be thrown away, and
recovered later by applying CNOT operations on a 0-state system.
In this way all of the information in the original quantum state, before performing the QSVD,
is packed into one quantum subsystem; plus the classical matrices Q and
V that are needed to specify the quantum transformations which would
recover the original state.
Continuing the example from above:
>> D
D =
0.94115 0 0 0
0 0.3076 0 0
0 0 0.10762 0
0 0 0 0.089686
>> resize(diag(D),4,4)
ans =
0.94115 0 0 0
0.3076 0 0 0
0.10762 0 0 0
0.089686 0 0 0
>> X=ans(:,1) # state vector for subsystem x
X =
0.94115
0.3076
0.10762
0.089686
>> Q'*diag(X,4,4)*V # equals A, the original bipartite (x,y) system
ans =
0.11312 0.22624 0.15082 0.30165
0.36198 0.13574 0.27148 0.18099
0.17237 0.34474 0.12928 0.25855
0.31027 0.20684 0.41369 0.15513
>> A
A =
0.11312 0.22624 0.15082 0.30165
0.36198 0.13574 0.27148 0.18099
0.17237 0.34474 0.12928 0.25855
0.31027 0.20684 0.41369 0.15513
In matrix form, the first CNOT moves the diagonal of D into
the first column, which then represents the state vector X
for subsystem x, and subsystem y is in the 0 state.
The other columns can then be thrown away. To recover the original system,
another CNOT puts X back on the diagonal, and then transformations
Q' and V undo the QSVD.
fidelity = (|<A,B>|/(||A||*||B||))2 where A (B) is the reduced state with qubit k equal to 0 (1), and <A,B>/(||A||*||B||)) is the cosine of the angle between vectors A and B.
If qubit k is not entangled with any of the other qubits, then the reduced state A should be the same as the reduced state B, i.e. the reduced state should be independent of qubit k. A and B represent the same state if the angle is 0 or 180 degrees, which corresponds to the cosine having a magnitude squared value of 1.
[1] Quantum Singular Value Decomposer, Carlos Bravo-Prieto, Diego Garcia-Martin, Jose I. Latorre, Phys. Rev. A 101, 062310 (2020), https://arxiv.org/abs/1905.01353
[2] Nash Singular Value Decomposition for Real or Complex, R. Perry, Jan 2022, https://fog.misty.com/perry/svd/svd.html
[3] No-cloning Theorem, https://en.wikipedia.org/wiki/No-cloning_theorem