Import["https://quest.qtechtheory.org/QuEST.m"]
env = CreateDownloadedQuESTEnv[];

(Subscript[H, 2]) Hamiltonian in BK representation, reduced

hamil = .2252 + .3435 (Subscript[Z, 0]) -.4347 (Subscript[Z, 1]) + .5716 (Subscript[Z, 0])(Subscript[Z, 1]) + .091 (Subscript[Y, 0])(Subscript[Y, 1]) + .091 (Subscript[X, 0])(Subscript[X, 1]);

ground = Min @ Eigenvalues @ CalcPauliSumMatrix @ hamil

A universal ansatz

numQb = 3;
numParams = 8;
ancilla = 2;
ansatz = Circuit[ 
	(Subscript[Ry, 1])[(Subscript[θ, 1])](Subscript[Ry, 0])[(Subscript[θ, 3])] (Subscript[Rz, 1])[(Subscript[θ, 2])] (Subscript[Rz, 0])[(Subscript[θ, 4])] (Subscript[C, 1])[(Subscript[X, 0])](Subscript[Ry, 1])[(Subscript[θ, 5])] (Subscript[Ry, 0])[(Subscript[θ, 7])] (Subscript[Rz, 1])[(Subscript[θ, 6])] (Subscript[Rz, 0])[(Subscript[θ, 8])]];

DrawCircuit[ansatz, numQb]

Li derivative circuits

insertGateDeriv[circ_, p_] := 
	circ /. {g:(Subscript[Ry, qb_])[(Subscript[θ, p])]:>{(Subscript[C, ancilla])[(Subscript[Y, qb])], g},g:(Subscript[Rz, qb_])[(Subscript[θ, p])]:>{(Subscript[C, ancilla])[(Subscript[Z, qb])],g}} // Flatten

getCircuitA[p1_,p2_] := 
	Join[{(Subscript[H, ancilla])}, insertGateDeriv[insertGateDeriv[ansatz,p1],p2], {(Subscript[H, ancilla])}]

getCircuitC[p_, h_] := 
	Join[
		{(Subscript[H, ancilla]), (Subscript[Rz, ancilla])[-π/2], (Subscript[X, ancilla])}, insertGateDeriv[ansatz,p], 
		{(Subscript[X, ancilla])},  If[h==1,{},(Subscript[C, ancilla])/@(List@@hamil[[h]])[[2;;]]], 
		{(Subscript[H, ancilla])}]

DrawCircuit[getCircuitA[3,5], numQb]

DrawCircuit[getCircuitC[4,3], numQb]

Manipulate[Column@{
	Row[{Style[ (BraKet[(("∂ψ")/("∂"(Subscript[θ, p1]))), (("∂ψ")/("∂"(Subscript[θ, p2])))]),Blue],": ", DrawCircuit[getCircuitA[p1,p2] ~Append~ (Subscript[M, ancilla]), numQb]}],
	Row[{Style[Row[{(Subscript["(Bra[ψ])H", h]),(Ket[(("∂ψ")/(Subscript["∂θ", p2]))])}],Red],": ", DrawCircuit[getCircuitC[p2,h] ~Append~ (Subscript[M, ancilla]), numQb]}]},
	{{p1,1,(Subscript[θ, i])},1,numParams,1}, {{p2,1,(Subscript[θ, j])},1,numParams,1}, {{h,1,(Subscript[H, j])},1,Length[hamil],1}, Paneled->False]

A circuits need not apply gates after entangling the ancilla

optimizeA[circ_] := 
	ReplaceRepeated[
		circ[[;;Position[circ,  (Subscript[C, ancilla])[_]][[-1,1]]]],
		{a___, g1:(Subscript[Except[C], qb1_])[_],g2:((Subscript[C, ancilla])[(Subscript[_, qb2_])]|PatternSequence[(Subscript[_, qb2_])[_], (Subscript[C, ancilla])[(Subscript[_, qb2_])]])} /; (qb1 ≠ qb2) :> {a,g2}
	] ~Append~ (Subscript[H, ancilla])

Manipulate[
	Row @ {Style[(BraKet[(("∂ψ")/("∂"(Subscript[θ, p1]))), (("∂ψ")/("∂"(Subscript[θ, p2])))]),Blue],": ",
		DrawCircuit[getCircuitA[p1,p2] ~Append~ (Subscript[M, ancilla]), numQb], Style[" -> ",Red],
		DrawCircuit[optimizeA@getCircuitA[p1,p2]~Append~ (Subscript[M, ancilla]), numQb]},
	{{p1,1,(Subscript[θ, i])},1,numParams,1}, {{p2,1,(Subscript[θ, j])},1,numParams,1}, Paneled->False]

Direct simulation of experimental routine

ψ = CreateQureg[3];
ϕ = CreateQureg[3];

getMatrixA[params_] := (1/4)Table[
	ApplyCircuit[ optimizeA@getCircuitA[r,c] /.params, InitZeroState[ψ]];
	2 CalcProbOfOutcome[ψ,ancilla,0]-1,
	{r,numParams}, {c,numParams}]

getVectorC[params_] := (1/2)Table[
	Sum[
		ApplyCircuit[getCircuitC[p,h]/.params, InitZeroState[ψ]];
		If[h==1,hamil[[h]],hamil[[h,1]]] 
		(2 CalcProbOfOutcome[ψ,ancilla,0]-1),
		{h,Length[hamil]}],
	{p,numParams}]

directImagTime[θi_, Δt_, steps_] := Module[
	{params, Δparams, matrA,vecC, energy, energies},
	params = θi;
	energies = {};
	Do[
		matrA = getMatrixA[params];
		vecC = getVectorC[params];
		Δparams = Δt LinearSolve[matrA, vecC, Method -> "Multifrontal"];
		params = MapIndexed[(Subscript[θ, #2[[1]]])->#1&,  params[[All,2]] + Δparams];

		ApplyCircuit[ansatz/.params, InitZeroState[ψ]];
		energy = CalcExpecPauliSum[ψ, hamil, ϕ];
		AppendTo[energies, energy];

		plotDirect = ListLinePlot[
			energies,
			GridLines-> {{},{{ground,Directive[Red,Dashed]}}},
			PlotRange-> {{0,steps},{1.05 ground,-ground}},
			AxesLabel-> {"Iterations", "Energy"}
		],
		steps
	]
]

Dynamic[plotDirect]
directImagTime[Table[(Subscript[θ, p])->RandomReal[{0,2π}],{p,numParams}], .1, 50]

Optimised simulation

dψ = CreateQuregs[numQb, numParams];
hψ = CreateQureg[numQb];

optimisedImagTime[θi_, Δt_, steps_] := Module[
	{params, Δparams, matrA,vecC, energy, energies},
	params=θi;
	energies = {};
	Do[
		CalcQuregDerivs[ansatz, InitZeroState[ψ], params, dψ];
		matrA = CalcInnerProducts[dψ] // Re;
	
		ApplyCircuit[ansatz/.params, InitZeroState[ψ]];
		ApplyPauliSum[ψ, hamil, hψ];
		vecC = -CalcInnerProducts[hψ, dψ] // Re;
	
		Δparams = Δt LinearSolve[matrA, vecC, Method -> "Multifrontal"];
		params = MapIndexed[(Subscript[θ, #2[[1]]])->#1&,  params[[All,2]] + Δparams];
	
		energy = CalcExpecPauliSum[ψ, hamil, ϕ];
		AppendTo[energies, energy],
		steps];
	energies]

optimisedRepeatImagTime[Δt_, steps_, reps_] :=Module[
	{evo, evos},
	evos = {};
	Do[
		evo = optimisedImagTime[Table[(Subscript[θ, p])->RandomReal[{0,2π}],{p,numParams}], Δt, steps];
		AppendTo[evos,evo];
		plotOptimised = ListLinePlot[
			evos,
			GridLines-> {{},{{ground,Directive[Red,Dashed]}}},
			PlotRange-> {{0,steps},{1.05 ground,-ground}},
			AxesLabel-> {"Iterations", "Energy"}
		],
		reps]]

Dynamic[plotOptimised]
optimisedRepeatImagTime[.05, 100, 10]