Harrow–Hassidim–Lloyd (HHL) algorithm

Install Wolfram quantum framework:
In[]:=
PacletInstall["https://www.wolfr.am/DevWQCF",ForceVersionInstall->True]​​<<Wolfram`QuantumFramework`
Out[]=
PacletObject
Name: Wolfram/QuantumFramework
Version: 1.6.2

Given a matrix A and a vector b, we want to solve
A·x=b
Generate HHL circuit that takes
A
,
b
, and number of ancillary qubits:
In[]:=
ClearAll[HHL]​​HHL[A_?SquareMatrixQ,b_?VectorQ,prec:_Integer?Positive:4]:=​​Block{d,qpe,t0,yTilde,λTilde,θλ,γ},​​d=Log2@Length[b];​​t0=Floor@πNorm[A];​​yTilde=If[0<=#<=
prec-1
2
,#,#-
prec
2
]&/@Range[0,
prec
2
-1];​​λTilde=
2π
t0
yTilde
prec
2
;​​γ=Min[DeleteCases[Abs[λTilde],0.|0]];​​θλ=If[#==0||#==0.,0.,2ArcSin[γ/#]]&/@λTilde;​​QuantumCircuitOperator[{​​QuantumState[Normalize[b]]->1+prec+Range[d],​​"Barrier",​​qpe=QuantumCircuitOperator["PhaseEstimation"[QuantumOperator[MatrixExp[t0IA]],prec],1+Range[prec+d]][[d+1;;-prec-1]],​​"Barrier",​​QuantumCircuitOperator["GrayOracle"[θλ,prec,1],"Eigenvalue Inversion"],​​"Barrier",​​qpe["Dagger"],​​"Barrier",​​QuantumOperator[QuantumState["1"<>StringJoin[ConstantArray["0",prec]]]["Dagger"]]​​}​​]/;IntegerQ[d]​​
Generate a complex-valued random Hermitian of
n
-qubits (a matrix of
n×n
)
In[]:=
n=4;​​a=QuantumOperator["RandomHermitian",n]["Matrix"];​​HermitianMatrixQ[a]
Out[]=
True
Generate a complex-valued random vector
b
:
In[]:=
b=QuantumOperator["RandomUnitary",n]["Matrix"][[All,1]];​​Norm[b]
Out[]=
1.
Create HHL circuit:
In[]:=
hhl=N@HHL[a,b,6];​​hhl["Diagram","ShowGateLabels"->False]
Out[]=
Given the result of above circuit, find corresponding amplitudes, and then normalize them to get solution of HHL solver:
In[]:=
xHHL=hhl[]["Normalize"]["AmplitudesList"]//Chop
Out[]=
{0.150198+0.266816,-0.306177-0.156923,0.177673+0.185801,-0.0630882-0.14418,-0.162317+0.24268,0.295957-0.105614,0.199304+0.0705461,0.0332224-0.0539288,0.324446+0.107075,0.00856717-0.106678,-0.189394+0.208001,-0.0599942-0.295949,-0.085392-0.0445195,-0.139345+0.193145,0.129758+0.117405,-0.167164+0.203132}
Find the corresponding solution using classical approach:
In[]:=
xClassical=Normalize@LinearSolve[a,b]
Out[]=
{0.147332+0.269409,-0.310906-0.155841,0.178197+0.191727,-0.0597219-0.134954,-0.169003+0.239835,0.301364-0.107247,0.194155+0.0756908,0.030282-0.0516326,0.332801+0.1021,0.00678392-0.0907113,-0.184115+0.202898,-0.0637374-0.290061,-0.0899423-0.0417513,-0.144571+0.192868,0.128466+0.118092,-0.16236+0.203951}
Compare classical result with the quantum one (using the inner product as a measure; the closer to one the better)"
In[]:=
Abs[Conjugate[xClassical].xHHL]^2
Out[]=
0.999141
As mentioned, the number of phase-estimation qubits (3rd argument) determines the precision of final outcome. Execute HHL with less phase-estimation qubits::
In[]:=
xHHLLowPrecision=N[HHL[a,b,5]][]["Normalize"]["AmplitudesList"]//Chop
Out[]=
{0.114006+0.325262,-0.323725-0.125372,0.15549+0.250044,-0.0920291+0.0429472,-0.295294+0.141047,0.346626-0.0697672,0.0923399+0.182321,-0.0481017-0.055415,0.430408-0.0650979,-0.0218042+0.161607,-0.0994849+0.186502,-0.0642308-0.0901146,-0.116693+0.00177034,-0.17157+0.140782,0.0745725+0.152716,-0.0225325+0.13835}
Quantify the good ness of HHL result in compared with the classical result:
In[]:=
Abs[Conjugate[xClassical].xHHLLowPrecision]^2
Out[]=
0.737561
Final comment: we put the initial
|b〉
state directly into our HHL circuit. We also have state-preparation circuit, meaning one can generate the state using some gates
In[]:=
prep=QuantumCircuitOperator[QuantumState[b],"StatePreparation"];​​prep["Diagram","ShowGateLabels"->False]
Out[]=
In[]:=
prep[]==QuantumState[b]
Out[]=
True