59 create equations of motion from a set of tensor indices (#60)

* fix notation for evolve

* add TensorProblem

* update struct

* add X tensors

* add tensor test

* update examples and docs

* fix notation

docs/src/examples/mean_field.md

@@ -48,6 +48,8 @@ m_i(t) = h_i + \sum_{j=1}^N J_{ij} n_j^z(t).
 \end{align}
 ```
 
+## Mean-Field Approximate Optimization
+
 To implement these dynamics within `QAOA.jl`, we begin by defining a schedule:
 ```julia
 using QAOA, LinearAlgebra
@@ -82,9 +84,9 @@ This is all we need to call the mean-field dynamics. The initial values are
 ```julia
 S = [[1., 0., 0.] for _ in 1:N-1]
 ```
-where we have taken into account that the _final spin is fixed_. The final vector of spins is then obtained as
+where we have taken into account that the _final spin is fixed_. The final vector of spins is then obtained via
 ```julia
-S = evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ)
+evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ)
 ```
 The energy expectation value in mean-field approximation is 
 ```julia
@@ -93,4 +95,56 @@ E = expectation(S, mf_problem.local_fields, mf_problem.couplings)
 and the solution of the algorithm can be retrieved by calling
 ```julia
 sol = mean_field_solution(S)
-```
+```
+
+## Equations of Motion via ODE Solver
+
+There is also the option to solve the full mean-field equations of motion (s. e.g. equation (9) [in this paper](https://arxiv.org/pdf/2403.11548)) via `OrdinaryDiffEq.jl`. For a simple linear annealing schedule, we can do this by calling
+```
+T_final = p * τ
+schedule_function = t -> t/T_final
+sol = evolve(mf_problem.local_fields, mf_problem.couplings, T_final, schedule_function)
+```
+
+### Tensor Problem Definition
+
+Instead of being restricted to simple QUBO problems defined by `local_fields` ``h_i`` and `couplings` ``J_{ij}`` as introduced above, we also want to be able to solve problems defined in terms of arbitrary tensors, e.g.
+```math
+\begin{align}
+    H(t) = (1 - s(t)) \sum_{i=1}^N n_i^x(t) + s(t)\sum_{i=1}^N \bigg[ J_i + \sum_{j>i} \big[ J_{ij} + \sum_{k>j} J_{ijk}n_k^z(t) + \cdots \big]  n_j^z(t) \bigg] n_i^z(t).
+\end{align}
+```
+To define the tensors for the standard QUBO set-up discussed here, we build the dictionaries
+```julia
+xtensor = Dict([(i, ) => 1.0 for i in 1:mf_problem.num_qubits])
+```
+for the local transverse-field driver and
+```julia
+ztensor = Dict()
+for (i, h_i) in enumerate(mf_problem.local_fields)
+    if h_i != 0.0
+        ztensor[(i,)] = h_i
+    end
+end
+
+for i in 1:mf_problem.num_qubits
+    for j in i+1:mf_problem.num_qubits
+        if mf_problem.couplings[i, j] != 0.0
+            ztensor[(i, j)] = mf_problem.couplings[i, j]
+        end
+    end
+end
+```
+for the problem Hamiltonian. We then feed these to
+```julia
+tensor_problem = TensorProblem(mf_problem.num_qubits, xtensor, ztensor)
+```
+and finally call
+```julia
+sol = evolve(tensor_problem, T_final, schedule_function)
+```
+which reproduces the same result as above.
+
+
+!!! note
+    `TensorProblem` can also deal with arbitrary higher-order tensors for the driver Hamiltonian!

docs/src/examples/prime_number.md

@@ -78,7 +78,7 @@ and call the mean-field AOA:
 mf_problem = Problem(p, -h, -J)
 S = [[1., 0., 0.] for _ in 1:N]
 
-S = evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ)
+evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ)
 ```
 As you can verify with this [notebook](https://github.com/FZJ-PGI-12/QAOA.jl/blob/master/notebooks/prime_number.ipynb), the solution returned by our algorithm agrees with the true solution:
 ```julia

docs/src/index.md

@@ -38,8 +38,8 @@ optimize_parameters(problem::Problem, beta_and_gamma::Vector{Float64}; niter::In
 ### Mean-Field Approximate Optimization Algorithm
 
 ```@docs
-evolve(S::Vector{<:Vector{<:Real}}, h::Vector{<:Real}, J::Matrix{<:Real}, β::Vector{<:Real}, γ::Vector{<:Real})
-evolve(S::Vector{<:Vector{<:Vector{<:Real}}}, h::Vector{<:Real}, J::Matrix{<:Real}, β::Vector{<:Real}, γ::Vector{<:Real})
+evolve!(S::Vector{<:Vector{<:Real}}, h::Vector{<:Real}, J::Matrix{<:Real}, β::Vector{<:Real}, γ::Vector{<:Real})
+evolve!(S::Vector{<:Vector{<:Vector{<:Real}}}, h::Vector{<:Real}, J::Matrix{<:Real}, β::Vector{<:Real}, γ::Vector{<:Real})
 evolve(h::Vector{<:Real}, J::Matrix{<:Real}, T_final::Float64, schedule::Function; rtol=1e-4, atol=1e-6)
 expectation(S::Vector{<:Vector{<:Real}}, h::Vector{<:Real}, J::Matrix{<:Real})
 mean_field_solution(problem::Problem, β::Vector{<:Real}, γ::Vector{<:Real})

notebooks/max_cut.ipynb

@@ -234,7 +234,7 @@
     "S = [[[1., 0., 0.] for _ in 1:N-1] for _ in 1:p+1]\n",
     "\n",
     "# evolution with history\n",
-    "S = evolve(S, mf_problem.local_fields, mf_problem.couplings, β, γ);"
+    "evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ);"
    ]
   },
   {

notebooks/mean_field.ipynb

@@ -101,7 +101,7 @@
     "S = [[1., 0., 0.] for _ in 1:N-1] # fix final spin (i.e. leave it out)\n",
     "\n",
     "# evolution\n",
-    "S = evolve(S, mf_problem.local_fields, mf_problem.couplings, β, γ)"
+    "evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ)"
    ]
   },
   {
@@ -123,7 +123,7 @@
     "S = [[[1., 0., 0.] for _ in 1:N-1] for _ in 1:p+1]\n",
     "\n",
     "# evolution with history\n",
-    "S = evolve(S, mf_problem.local_fields, mf_problem.couplings, β, γ);"
+    "evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ);"
    ]
   },
   {
@@ -195,26 +195,96 @@
     "mean_field_solution(S[end])"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "a9beaf4c",
+   "metadata": {},
+   "source": [
+    "## Using the ODE solver"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a9beaf4c",
+   "id": "b9551394",
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "T_final = p*τ\n",
+    "schedule_function = t -> t/T_final\n",
+    "sol = QAOA.evolve(mf_problem.local_fields, mf_problem.couplings, T_final, schedule_function)\n",
+    "sol.u[end]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1ac23ab4",
+   "metadata": {},
+   "source": [
+    "### Using the tensor problem definition\n",
+    "\n",
+    "Note that `TensorProblem` can also deal with arbitrary higher-order tensors."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fe091fc9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "    # create tensors\n",
+    "    xtensor = Dict([(i, ) => 1.0 for i in 1:mf_problem.num_qubits])\n",
+    "    ztensor = Dict()\n",
+    "    for (i, h_i) in enumerate(mf_problem.local_fields)\n",
+    "        if h_i != 0.0\n",
+    "            ztensor[(i,)] = h_i\n",
+    "        end\n",
+    "    end\n",
+    "    \n",
+    "    for i in 1:mf_problem.num_qubits\n",
+    "        for j in i+1:mf_problem.num_qubits\n",
+    "            if mf_problem.couplings[i, j] != 0.0\n",
+    "                ztensor[(i, j)] = mf_problem.couplings[i, j]\n",
+    "            end\n",
+    "        end\n",
+    "    end "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2e33b302",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define tensor problem\n",
+    "tensor_problem = TensorProblem(mf_problem.num_qubits, xtensor, ztensor)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15aa5b7d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol = QAOA.evolve(tensor_problem, T_final, schedule_function)\n",
+    "sol.u[end]"
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.8.2",
+   "display_name": "Julia 1.9.4",
    "language": "julia",
-   "name": "julia-1.8"
+   "name": "julia-1.9"
   },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.8.2"
+   "version": "1.9.4"
   }
  },
  "nbformat": 4,

notebooks/prime_number.ipynb

@@ -136,30 +136,9 @@
    "execution_count": 6,
    "id": "def6674f",
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "10-element Vector{Vector{Float64}}:\n",
-       " [0.32825975222477094, 0.11978494573299431, 0.9369616330699301]\n",
-       " [0.70620595460824, 0.14988113816534768, -0.6919601101928738]\n",
-       " [0.15180019423091565, -0.04152395891726101, 0.9875385875333125]\n",
-       " [0.2271263818640098, 0.011061686135514353, 0.9738024675267594]\n",
-       " [-0.32149023287393297, -0.07047142076523585, 0.9442869315107548]\n",
-       " [0.647238279110526, -0.7140992021808273, -0.26673008735187687]\n",
-       " [0.14563391972076617, -0.2647005386966897, 0.9532703636641859]\n",
-       " [-0.03218817700881943, -0.12424388710135918, 0.9917294882067359]\n",
-       " [-0.09497156319945207, 0.05694668862360597, -0.99384982609963]\n",
-       " [0.6105715312011801, 0.12405980734829188, 0.7821838463477425]"
-      ]
-     },
-     "execution_count": 6,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "S = evolve(S, mf_problem.local_fields, mf_problem.couplings, β, γ)"
+    "evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ)"
    ]
   },
   {
@@ -188,15 +167,15 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.8.2",
+   "display_name": "Julia 1.9.4",
    "language": "julia",
-   "name": "julia-1.8"
+   "name": "julia-1.9"
   },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.8.2"
+   "version": "1.9.4"
   }
  },
  "nbformat": 4,

notebooks/sherrington_kirkpatrick.ipynb

@@ -199,7 +199,7 @@
     "S = [[[1., 0., 0.] for _ in 1:N-1] for _ in 1:p+1]\n",
     "\n",
     "# evolution with history\n",
-    "S = evolve(S, mf_problem.local_fields, mf_problem.couplings, β, γ);"
+    "evolve!(S, mf_problem.local_fields, mf_problem.couplings, β, γ);"
    ]
   },
   {

src/QAOA.jl

@@ -10,7 +10,7 @@ using OrdinaryDiffEq
 import Distributions, Random
 
 include("problem.jl")
-export Problem
+export Problem, TensorProblem
 
 include("optimization.jl")
 export cost_function, optimize_parameters
@@ -21,7 +21,7 @@ export sherrington_kirkpatrick, partition_problem, max_cut, min_vertex_cover
 include("circuit.jl")
 
 include("mean_field.jl")
-export evolve, expectation, mean_field_solution
+export evolve!, evolve, expectation, mean_field_solution
 
 include("fluctuations.jl")
 export evolve_fluctuations

src/fluctuations.jl

@@ -56,7 +56,7 @@ function evolve_fluctuations(problem::Problem, τ::Real, β::Vector{<:Real}, γ:
     # evolution
     S = [[1., 0., 0.] for _ in 1:num_qubits]
     S = [S for _ in 1:num_layers+1]
-    S = evolve(S, local_fields, couplings, β, γ)
+    evolve!(S, local_fields, couplings, β, γ)
     solutions = mean_field_solution(S[end])  
 
     M = 1.0I(2num_qubits)