-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathAdiabaticTimeEvolution.m
336 lines (299 loc) · 14 KB
/
AdiabaticTimeEvolution.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
% ------------------------------------------------------------------------------
% Adiabatic Time Evolution for the TFXY Hamiltonian. The parameters of this
% script are setup to compute the circuits used for the adiabatic state
% preparation experiment presented in Fig. 6 of:
% Algebraic Compression of Quantum Circuits for Hamiltonian Evolution,
% Efekan Kockcu, Daan Camps, Lindsay Bassman, J. K. Freericks,
% Wibe A. de Jong, Roel Van Beeumen, and Alexander F. Kemper (2021)
% ------------------------------------------------------------------------------
% ### Parameters to adjust #####################################################
% Overall parameters -----------------------------------------------------------
N = 5; % number of spins in the TFXY Hamiltonian
nsteps = 240; % number of Trotter steps:
% * 240 for simulation with large Trotter
% error (purple line in Fig. 6)
% * 1200 for simulation with small Trotter
% error (red line in Fig. 6)
dt = 0.25; % size of a Trotter step:
% * 0.25 for simulation with large Trotter
% error (purple line in Fig. 6)
% * 0.05 for simulation with small Trotter
% error (purple line in Fig. 6)
% circuit storage directory
folder = 'ADD LOCAL PATH HERE';
saveCircuit = false; % true/false: save quantum circuits to file.
% ------------------------------------------------------------------------------
% Transverse field hz ----------------------------------------------------------
hz.start = -1.0; % value of hz at start of simulation
hz.end = -1.0; % value of hz at end of simulation
hz.startRamp = 1; % step when ramping starts
hz.endRamp = nsteps; % step when ramping ends
hz.rampType = @linearRamp; % type of ramp: linear or sigmoid
% ------------------------------------------------------------------------------
% Coupling parameter Jx --------------------------------------------------------
Jx.start = 0.0; % value of Jx at start of simulation
Jx.end = -2.0; % value of Jx at end of simulation
Jx.startRamp = 1; % step when ramping starts
Jx.endRamp = round( 0.5 * nsteps ); % step when ramping ends
Jx.rampType = @linearRamp; % type of ramp: linear or sigmoid
% ------------------------------------------------------------------------------
% Coupling parameter Jy --------------------------------------------------------
Jy.start = 0.0; % value of Jy at start of simulation
Jy.end = 0.0; % value of Jy at end of simulation
Jy.startRamp = 1; % step when ramping starts
Jy.endRamp = nsteps; % step when ramping ends
Jy.rampType = @linearRamp; % type of ramp: linear or sigmoid
% ------------------------------------------------------------------------------
% ##############################################################################
clc
fprintf( 1, '\n---------------------------------------------------------------\n');
fprintf( 1, 'Generating and validating quantum circuits for TFXY Hamiltonian\n' );
fprintf( 1, '---------------------------------------------------------------\n\n');
% generate all parameter arrays for full simulation.
hz_arr = arrayFromCell( hz, nsteps );
Jx_arr = arrayFromCell( Jx, nsteps );
Jy_arr = arrayFromCell( Jy, nsteps );
% set timing arrays.
tfinal = ( nsteps - 1 ) * dt;
t_arr = linspace( 0, tfinal, nsteps );
% print information about simulation to command window.
fprintf( 1, 'Number of spins = \t\t\t%d\n', N );
fprintf( 1, 'Number of Trotter steps = \t\t%d\n', nsteps );
fprintf( 1, 'Stepsize of Trotter step = \t\t%.3f\n', dt );
fprintf( 1, 'Final simulation time = \t\t%.3f\n', tfinal );
% save some plot colors.
color1 = [252, 186, 3]./255;
color2 = [3, 92, 59]./255;
color3 = [167, 38, 222]./255;
color4 = [37, 122, 118]./255;
% plot the simulation parameter curves.
figure(1); clf
plot(t_arr, Jx_arr, 'color', color1, 'LineWidth', 2, 'DisplayName', 'J_x');
hold on
plot(t_arr, Jy_arr, 'color', color2, 'LineWidth', 2,'DisplayName', 'J_y');
plot(t_arr, hz_arr, 'color', color3, 'LineWidth', 2,'DisplayName', 'h_z');
xlabel('Simulation time');
title('Simulation parameters for TFXY');
legend();
% compute the Hamiltonian eigenvalues throughout simulation.
energies = zeros(nsteps, 2^N);
for i = 1:nsteps
energies(i, :) = eig( hamTFXY(N, Jx_arr(i), Jy_arr(i), hz_arr(i)) );
end
% compute the magnetization observable O
Z = qclab.qgates.PauliZ.matrix;
O = zeros(2^N, 2^N);
for i = 0:N-1
O = O + kron(kron( qclab.qId(i), Z), qclab.qId(N - i - 1) );
end
% start state(s) for simulation.
psi0 = zeros(2^N, 1);
psi0(end) = 1; % all 1's initial state.
psi_ex = psi0;
psi_tr = psi0;
% setup arrays for results.
% energy
energy_ex = zeros( nsteps, 1 ); % energy with exact unitary
energy_tr = zeros( nsteps, 1 ); % energy with trotter unitary
energy_fc = zeros( nsteps, 1 ); % energy with full circuit
energy_cc = zeros( nsteps, 1 ); % energy with compressed circuit
% magnetization
mag_ex = zeros( nsteps, 1 ); % magnetization with exact unitary
mag_tr = zeros( nsteps, 1 ); % magnetization with trotter unitary
mag_fc = zeros( nsteps, 1 ); % magnetization with full circuit
mag_cc = zeros( nsteps, 1 ); % magnetization with compressed circuit
% do the full simulation and generate + save circuits.
fullCircuitMat = qclab.QCircuit( N );
fullCircuitRot = qclab.QCircuit( N );
compressedCircuitRot = qclab.QCircuit( N );
TFXY = @f3c.qgates.RotationTFXY;
numerical_error = zeros( nsteps, 1 ); % numerical error on compression algorithm
circLayers = 1;
for i = 1:nsteps
% Print info to command window.
fprintf( 1, 'Current time = \t\t\t%.3f\t(%d/%d)\n', t_arr(i), i, nsteps );
% Hamiltonian at this timestep.
H = hamTFXY( N, Jx_arr(i), Jy_arr(i), hz_arr(i) );
% full unitaries: exact and trotter ------------------------------------------
psi_ex = expm( -1i * dt * H ) * psi_ex ;
psi_tr = expm( -1i * dt * hamTFXY( N, Jx_arr(i), Jy_arr(i), 0 ) ) * ...
expm( -1i * dt * hamTFXY( N, 0, 0, hz_arr(i) ) ) * psi_tr ;
% energies
energy_ex(i) = real( psi_ex' * H * psi_ex );
energy_tr(i) = real( psi_tr' * H * psi_tr );
% magnetization
mag_ex(i) = -real( psi_ex' * O * psi_ex )/N;
mag_tr(i) = -real( psi_tr' * O * psi_tr )/N;
% quantum circuits -----------------------------------------------------------
[timeStepMat, timeStepRot] = f3c.timeStepTFXY( N, dt, Jx_arr(i), Jy_arr(i), hz_arr(i) );
if circLayers < ceil(N/2) % not yet at minimal depth size
for j = 1:timeStepMat.nbGates
fullCircuitMat.push_back( timeStepMat.gates(j) );
fullCircuitRot.push_back( timeStepRot.gates(j) );
end
squareCircuit = copy( fullCircuitMat );
elseif circLayers == ceil(N/2) % becomes minimal depth at current size
if mod( N, 2 ) == 0 % even number of spins
for j = 1:timeStepMat.nbGates
fullCircuitMat.push_back( timeStepMat.gates(j) );
fullCircuitRot.push_back( timeStepRot.gates(j) );
end
squareCircuit = f3c.SquareCircuit( fullCircuitMat );
triangleCircuit = f3c.TriangleCircuit( squareCircuit );
else
% add first vertical layer to make it square
for j = 1:(N-1)/2
fullCircuitMat.push_back( timeStepMat.gates(j) );
fullCircuitRot.push_back( timeStepRot.gates(j) );
end
squareCircuit = f3c.SquareCircuit( fullCircuitMat );
triangleCircuit = f3c.TriangleCircuit( squareCircuit );
% add second vertical layer to fullCircuit and merge with triangle
for j = ((N-1)/2)+1:timeStepMat.nbGates
fullCircuitMat.push_back( timeStepMat.gates(j) );
triangleCircuit.merge( timeStepMat.gates(j), 'R' );
fullCircuitRot.push_back( timeStepRot.gates(j) );
end
squareCircuit = f3c.SquareCircuit( triangleCircuit );
end
else % larger than minimal depth
for j = 1:timeStepMat.nbGates
fullCircuitMat.push_back( timeStepMat.gates(j) );
triangleCircuit.merge( timeStepMat.gates(j), 'R' );
fullCircuitRot.push_back( timeStepRot.gates(j) );
end
squareCircuit = f3c.SquareCircuit( triangleCircuit );
end
% rotation representation of compressed circuit.
compressedCircuitRot.clear();
for j = 1:squareCircuit.nbGates
compressedCircuitRot.push_back( TFXY( squareCircuit.gates(j) ) );
end
% compute the state based on the quantum circuits
psi_fc = fullCircuitMat.apply('R', 'N', N, psi0 );
psi_cc = squareCircuit.apply('R', 'N', N, psi0 );
% energies
energy_fc(i) = real( psi_fc' * H * psi_fc );
energy_cc(i) = real( psi_cc' * H * psi_cc );
% magnetization
mag_fc(i) = -real( psi_fc' * O * psi_fc )/N;
mag_cc(i) = -real( psi_cc' * O * psi_cc )/N;
numerical_error(i) = ...
norm( fullCircuitRot.matrix - compressedCircuitRot.matrix, 'fro' );
if mod(i, 10) == 0
fprintf(1, 'Numerical error on compressed circuit: %e\n', ...
numerical_error(i) );
end
% save the circuits to file.
if saveCircuit
fname_fc = [folder,'TFXY_full_circuit_t_',num2str(t_arr(i)),'.qasm'];
fid = fopen( fname_fc, 'w+' );
fprintf( fid, '// Generated by f3c\n');
fprintf( fid, '// https://github.com/QuantumComputingLab/f3c\n\n' );
fprintf( fid, '// Current time = \t%6.4f [dt = %6.4f][timestep = %d/%d]\n', ...
t_arr(i), dt, i, nsteps );
fprintf( fid, '// Jx = %6.4f,\t Jy = %6.4f,\t hz = %6.4f\n\n', ...
Jx_arr(i), Jy_arr(i), hz_arr(i) );
fprintf( fid, 'OPENQASM 2.0;\ninclude "qelib1.inc";\n\n');
fprintf( fid, 'qreg q[%d];\n',N);
fullCircuitRot.toQASM( fid );
fclose( fid );
fname_cc = [folder,'TFXY_compressed_circuit_t_',num2str(t_arr(i)),'.qasm'];
fid = fopen( fname_cc, 'w+' );
fprintf( fid, '// Generated by f3c\n');
fprintf( fid, '// https://github.com/QuantumComputingLab/f3c\n\n' );
fprintf( fid, '// Current time = \t%6.4f [dt = %6.4f][timestep = %d/%d]\n', ...
t_arr(i), dt, i, nsteps );
fprintf( fid, '// Jx = %6.4f,\t Jy = %6.4f,\t hz = %6.4f\n\n', ...
Jx_arr(i), Jy_arr(i), hz_arr(i) );
fprintf( fid, 'OPENQASM 2.0;\ninclude "qelib1.inc";\n\n');
fprintf( fid, 'qreg q[%d];\n',N);
compressedCircuitRot.toQASM( fid );
fclose( fid );
end
circLayers = circLayers + 1;
end
% plot the Hamiltonian eigenvalues throughout simulation together
% with approximate ground states.
figure(2); clf
axEnergies = axes;
p = plot(t_arr, energies, 'k-');
for i = 1:length(p)
p(i).Annotation.LegendInformation.IconDisplayStyle = 'off';
end
hold on
plot( t_arr, energy_ex, 'color', color1, 'LineWidth', 2, 'DisplayName', 'Exact Unitary' );
plot( t_arr, energy_tr, 'color', color2, 'LineWidth', 2, 'DisplayName', 'Trotter Unitary');
plot( t_arr, energy_fc, '--','color', color3, 'LineWidth', 2, 'DisplayName', 'Full Circuit');
plot( t_arr, energy_cc, ':','color', color4, 'LineWidth', 2, 'DisplayName', 'Compressed Circuit');
xlabel('Simulation time');
title('Energies during TFXY simulation');
legend();
% plot the magnetization.
figure(3); clf
plot( t_arr, mag_ex, 'color', color1, 'LineWidth', 2, 'DisplayName', 'Exact Unitary' );
hold on
plot( t_arr, mag_tr, 'color', color2, 'LineWidth', 2, 'DisplayName', 'Trotter Unitary');
plot( t_arr, mag_fc, '--','color', color3, 'LineWidth', 2, 'DisplayName', 'Full Circuit');
plot( t_arr, mag_cc, ':','color', color4, 'LineWidth', 2, 'DisplayName', 'Compressed Circuit');
xlabel('Simulation time');
title('Magnetization during TFXY simulation');
legend();
% plot the numerical error on the compression algorithm
figure(4); clf
loglog( t_arr, numerical_error, 'color', color1, 'LineWidth', 2 );
xlabel('Simulation time');
title('Numerical error on compression algorithm');
ylabel('Frobenius norm of error');
ylim([1e-16,1e0])
% -- Functions ---------------------------------------------------------------
function f = arrayFromCell( in, nsteps )
f = in.rampType( in.startRamp, in.endRamp, nsteps );
f = in.start + f .* ( in.end - in.start );
end
% generates a linear ramp from 0 to 1 with npoints.
% it starts ramping up from point start \in [1, npoints-1]
% and ends ramping up at point end \in [2, npoints]
function f = linearRamp( startpoint, endpoint, npoints )
f = zeros( npoints, 1 );
ramp = linspace( 0, 1, endpoint - startpoint + 1);
f(startpoint:endpoint) = ramp;
f(endpoint+1:end) = 1;
end
% generates a sigmoid ramp from 0 to 1 with npoints
% it starts ramping up from point start \in [1, npoints-1]
% and ends ramping up at point end \in [2, npoints]
function f = sigmoidRamp( startpoint, endpoint, npoints )
f = zeros( npoints, 1 );
x = linspace( -10, 10, endpoint - startpoint + 1 );
f(startpoint:endpoint) = 1./( 1 + exp(-x) );
f(endpoint+1:end) = 1;
end
% generate the TFXY Hamiltonian for ordered (scalar) and disordered models
% (arrays)
function H = hamTFXY(N, Jx, Jy, hz)
assert( qclab.isNonNegInteger( N ) );
if isscalar( Jx )
Jx = repmat( Jx, N-1, 1 );
end
assert( length(Jx) == N-1 );
if isscalar( Jy )
Jy = repmat( Jy, N-1, 1 );
end
assert( length(Jy) == N-1 );
if isscalar( hz )
hz = repmat( hz, N, 1 );
end
assert( length(hz) == N );
XX = kron( qclab.qgates.PauliX.matrix, qclab.qgates.PauliX.matrix );
YY = kron( qclab.qgates.PauliY.matrix, qclab.qgates.PauliY.matrix );
Z = qclab.qgates.PauliZ.matrix;
H = zeros(2^N,2^N);
for i=0:N-2
H = H + Jx(i+1) * kron(kron( qclab.qId(i), XX ), qclab.qId(N - i - 2) ) ...
+ Jy(i+1) * kron(kron( qclab.qId(i), YY ), qclab.qId(N - i - 2) );
end
for i=0:N-1
H = H + hz(i+1) * kron(kron( qclab.qId(i), Z), qclab.qId(N - i - 1) );
end
end