Symbolic Quantum Optimization: 
Exploring QAOA-in-QAOA for Max-Cut Problem

Sebastián Rodríguez-Falcón
Wolfram Research South America

#WolframTechConf

ABSTRACT

Introduction to the Quantum Approximate Optimization Algorithm (QAOA) and QUBO formulation through the Max-Cut problem, using the Wolfram Quantum Framework for fully symbolic and detailed explanations accessible to audiences without a physics background, including the advanced QAOA-in-QAOA method.

Max-Cut problem

The Max–Cut problem is a combinatorial optimization problem that consists of finding a partition of the vertices of a graph into two complementary subsets in such a way that the total weight of the edges crossed by the cut is maximized.

Grid[{{Text@Classical Formulation,Text@Quantum Formulation},{HoldForm[],},{Style[C=(1/2)(∑_(𝒾,𝒿 ∈ V)) (1-(Subscript[s, 𝒾])(Subscript[s, 𝒿])), 
s∈{-1,+1},FontSize->20],Style[(Subscript[H, C])=(1/2)(∑_(𝒾,𝒿 ∈ V))(𝟙-(Subscript[Z, 𝒾])(Subscript[Z, 𝒿]))  
(Subscript[Z, i])(Subscript[Z, j])(Ket[(Subscript[x, 0])(Subscript[…x, n])])=(((-1))^(Subscript[x, i]))(((-1))^(Subscript[x, j]))(Ket[(Subscript[x, 0])(Subscript[…x, n])]),(Subscript[x, i])∈{0,1} ,FontSize->20]},{,}},Frame->All]

Quantum Approximate Optimization Algorithm

The Quantum Approximate Optimization Algorithm (QAOA), introduced by E. Farhi as a Variational Quantum Algorithm for the Max-Cut problem, alternates between cost and mixer unitaries to variationally minimize a cost Hamiltonian that encodes the objective function of a combinatorial optimization problem.

The QAOA process involves the following steps:

Defining a Cost Hamiltonian (Subscript[H, C]) and Mixer Hamiltonian (Subscript[H, M])

Defining the unitary gates: (Subscript[U, C])(γ)=(e^(-i*γ*(Subscript[H, C]))), (Subscript[U, M])(θ)=(e^(-i*θ*(Subscript[H, M])))

Applying the gates alternately in layers, forming the unitary: U(γ,θ)=(∏_(i=1)^N)(Subscript[U, c])(γ)(Subscript[U, M])(θ)

Preparing an Initial State and Applying U(t,θ)

Optimizing Parameters and Measuring

Understanding QAOA

Load Wolfram Quantum Framework

<<Wolfram`QuantumFramework`

Load Quantum Optimization Paclet

<<Wolfram`QuantumFramework`QuantumOptimization`

We are looking for “periodic” bitstrings:

miniQAOA=;
miniQAOA["Diagram",]

Identifying a periodic and regular bitstring:

superposition=miniQAOA[[1]];
Row[{#["Diagram"],#["Formula"]}]&@superposition

qc=miniQAOA[[;;2]];Row[{#["Diagram"],#["Formula"]}]&@qc

qc["Amplitudes"]

Plot[Evaluate[Im[Values[#1]]],{γ,-2 π,2 π},]&@qc["Amplitudes"][[;;2]]

qc["Probability"]

Plot[Evaluate[Values[#1]/. Im[γ]:>γ],{γ,-2 π,2 π},]&@qc["Probability"][[;;2]]

Mixing for every qubit:

miniQAOA["Diagram",]

miniQAOA["Simplify"]["Table"]

Plot3D[Evaluate[ComplexExpand@Values[#]],{γ,-2π,2π},{θ,-π,π},]&@miniQAOA[]["Probabilities"][[;;2]]

Example: 6-Node Graph

Find the Max-Cut for the following graph:

g6=Graph[]

In this more challenging example, we have a graph with six nodes, which results in a large cost Hamiltonian:

index6=MapApply[List,EdgeRules[g6]]/.x_Integer:>{x,x};
op𝟙6=StringRepeat["I",VertexCount[g6]];
sum6=Total[(1/2)(op𝟙6-StringReplacePart[op𝟙6,"Z",#])&/@index6]

Hc6=QuantumOperator[sum6]

Implement directly the QAOA circuit using the Graph functionalities:

qaoa6=QuantumCircuitOperator[{
"0"->VertexList[g6],
"H"->VertexList[g6],
Sequence@@Map[{"R",γ,"ZZ"}->#&,List@@@EdgeList[g6]],
"RX"[θ]->VertexList[g6]
},
"Label"->None,"Parameters"->{θ,γ}];

qaoa6["Diagram"]

Implement the symbolic cost function using the Cost Hamiltonian and the resulting parametrized state from the previous circuit:

(Ket[ϕ])=qaoa6[];(Bra[ϕ])=(SuperDagger[(Ket[ϕ])]);

costfunction6="Scalar"//(Bra[ϕ])@Hc6@(Ket[ϕ])//ComplexExpand//Simplify

Use NMaximize for the optimization:

result6=NMaximize[costfunction6,{θ,γ}]

We obtain multiple possible solutions; however, with just a single QAOA layer, we observe four solutions with higher probability:

(qaoa6@@ Values[Last[result6]])["ProbabilityPlot", ]

Solve it clasically:

FindMaximumCut[g6]

We can analyze the four most probable quantum solutions indicated by the ProbabilityPlot:

ReverseSortBy[Normal[(qaoa6 @@ Values[Last[result6]])["Probabilities"]],#[[2]]&][[;;4]]

Following the labels in the vertex, the solution (Ket[010011]) is the one found by the classical method. The following diagram corresponds to this solution:

maxcut=HighlightGraph[g,(Subgraph[g,#1]&)/@Last[FindMaximumCut[g]],VertexLabels->];cut=Show[{maxcut,},];
GraphicsRow[{g,cut},ImageSize->Large]

QAOA-in-QAOA

𝓍1=UndirectedEdge@@@{{1,2},{2,3},{3,4},{2,4},{1,3}};
𝓍2=UndirectedEdge@@@(4+{{1,2},{2,3},{3,4},{1,4},{2,4},{1,3}});
𝓍3=UndirectedEdge@@@(8+{{1,2},{2,3},{3,1}});
𝓍4=UndirectedEdge@@@(11+{{1,2}});
connection=UndirectedEdge@@@{{4,11},{10,7},{9,12},{13,5},{3,7},{1,8}};
edges={𝓍1,𝓍2,𝓍3,𝓍4};
fullEdges=Flatten[{{𝓍4,𝓍2,𝓍1,𝓍3},connection}];
subgraphs=Graph[#,EdgeWeight->Thread[Rule[#,1]],]&/@edges;
gsub=HighlightGraph[fullEdges,subgraphs,];
g13=Graph[fullEdges,PlotTheme->"NameLabeled",EdgeWeight->Thread[fullEdges->1]];
GraphicsRow[{g13,gsub},ImageSize->Full]

Divide Phase: Subgraphs

subweights=<|Thread[List@@@EdgeList[#]->PropertyValue[#,EdgeWeight]]|>&/@subgraphs;
subHc=Total[Table[(1/2)#[i]*(QuantumOperator["I",i]-QuantumOperator["Z",i]),{i,Keys@#}]]&/@subweights

Symbolic Variational Quantum Circuits

subQAOA=&/@subgraphs;

GraphicsGrid[Partition[#["Diagram","ShowEmptyWires"->False]&/@subQAOA,2],ImageSize->Full]

Cost Functions

QAOAstates=QuantumOperator[#]["State"]&/@subQAOA;

subcost=Module[{ket,bra,hc},;(bra.hc.ket)//ComplexExpand//FullSimplify]&/@Thread[{subHc,QAOAstates}];

contours=ContourPlot[#,{θ,-π,π},{γ,-π,π},]&/@subcost;

sg=Subgraph[gsub,#,ImageSize->Small]&/@{Range[4],Range[5,8],{9,10,11},{12,13}};
minigrids=Grid[{{#[[1]],#[[2]]},{Text[#[[3]],FormatType->TraditionalForm],SpanFromLeft}},Frame->All,Background->{White,{None,Lighter[Yellow,.9]}},ItemStyle->#[[4]]]&/@Thread[{sg,contours,subcost,{8,12,12,12}}];
Grid[{#[[{1,2}]],#[[{3,4}]]}&@minigrids,Frame->All,BaselinePosition->Baseline,Alignment->Top,Background->Lighter[Blend[{Red,Lighter[Black,0.7]}],.3]]

Optimization

subres=NMaximize[#,{θ,γ}]&/@subcost;

result13=(#1[[1]]@@Values[Last[#1[[2]]]]&)/@Thread[{subQAOA,subres}];
GraphicsGrid[Partition[#["ProbabilityPlot", ]&/@result13,2]]

states=Keys@ReverseSortBy[#["Probabilities"],#[[2]]]&/@result13;
groups=Map[states[[Sequence@@##1]][[1]]&,{{1,2},{2,6},{3,1},{4,1}}];
solution=Flatten[Thread/@Thread[(First[#1["Order"]]&)/@subQAOA->groups]];
minig=Graph[g13,VertexLabels->MapAt[Placed[#1,Center]&,solution,{All,2}]];incompleteHigh=HighlightGraph[minig,(Subgraph[g13,#1]&)/@GatherBy[solution,Last][[All,All,1]],VertexLabelStyle->Directive[Black,12],GraphHighlightStyle->"Thick"];
GraphicsRow[{gsub,incompleteHigh},ImageSize->Full]

Conquer Phase: Big Graph

In the Conquer phase, we construct a new graph where each node represents a subgraph, and the edges are based on their prior connections:

as=<|solution|>;
bigEdges=Complement[EdgeList[g13],Flatten[EdgeList/@Map[Subgraph[g13,#]&,edges]]];
smallToBig=Flatten[Thread/@Thread[VertexList/@{𝓍1,𝓍2,𝓍3,𝓍4}->Range[4]]];
bigedges=EdgeList[bigEdges]/.z:UndirectedEdge[x_,y_]:>Rule[z/.smallToBig,(((-1))^(as[x]))*(((-1))^(as[y]))]//Merge[#,Total]&//Normal;
g2=Graph[Reverse@#[[All,1]],EdgeWeight->#,EdgeLabels->"EdgeWeight",EdgeLabelStyle->20,PlotTheme->"NameLabeled",VertexStyle->{1 -> Hue[0,0.57,0.78,0.4], 2 ->Hue[0.15,0.82,0.87,0.4],3->Hue[0.75,0.5,0.66,0.4],4->Hue[0.05,0.62,0.87,0.6]},ImageSize->Medium]&@bigedges;
ggg=Show[{Graphics[{
Hue[0.15,0.82,0.87,0.4],Disk[{1.7,0.5},0.8],
Hue[0,0.57,0.78,0.4],Disk[{0.45,1.8},{0.8,1}],
Hue[0.75,0.5,0.66,0.4],Disk[{2.43,2.5},0.75],
Hue[0.05,0.62,0.87,0.6],Disk[{3.65,1.4},{0.3,1}]
}],incompleteHigh
}];
GraphicsRow[{ggg,g2},ImageSize->(10^2.9)]

We then assign weights based on inter-subgraph interactions and repeat the QAOA process at this higher level:

subweights2=<|Thread[List@@@EdgeList[#]->((1-PropertyValue[#,EdgeWeight])/2)]|>&@g2;

subHc2=Total[Table[(1/2)#[i]*(QuantumOperator["I",i]-QuantumOperator["Z",i]),{i,Keys@#}]]&@subweights2;

subQAOA2=QuantumCircuitOperator[{"+"->VertexList[#1],Sequence@@({"R",γ,"ZZ"}->#1&)/@List@@@EdgeList[#1],"RX"[θ]->VertexList[#1]},
"Label"->None,"Parameters"->{θ,γ}]&@g2;
subQAOA2["Diagram"]

QAOAstates2=QuantumOperator[#]["State"]&@subQAOA2;
subcost2=("Scalar"//(SuperDagger[QAOAstates2])@subHc2["State"]@QAOAstates2)//ComplexExpand//Simplify

subres2=NMaximize[subcost2,{θ,γ}]

Finally, we observe the resulting probability distribution and solve the MaxCut problem for the higher-level graph:

result2=subQAOA2["State"][<| Last[subres2]|>];

probplot=result2["ProbabilityPlot", ];

groups2=Part[Keys@ReverseSortBy[#["Probabilities"],#[[2]]],2][[1]]&@result2;
solution2=Flatten[Thread/@Thread[First[#["Order"]]&@subQAOA2->groups2]];
grouped2=HighlightGraph[g2,Map[Subgraph[g13,#]&,GatherBy[solution2,Last][[All,All,1]]],];
GraphicsRow[{probplot,grouped2},ImageSize->Large]

In the final step, we take the subgraphs that are yellow nodes in the higher level graph and flip their node values. For instance, in some subgraphs, yellow nodes may become red, as the global optimization balances local and global constraints.

flip=Keys[Select[solution2,MatchQ[#,Rule[_,1]]&]];
newsolution=#->1-as[#]&/@Keys[Select[smallToBig,MatchQ[#,Rule[_,Alternatives@@flip]]&]];
finalsolution=Normal[Association[Join@@{solution,newsolution}]];
finalsol=HighlightGraph[g13,Map[Subgraph[g13,#]&,GatherBy[finalsolution,Last][[All,All,1]]],];
GraphicsRow[{incompleteHigh,Style["->",20],finalsol},ImageSize->Full]

Thanks!

Sebastián Rodríguez-Falcón

References

“Quantum Optimization Tech Note”. Rodriguez, S. (2024). Wolfram Language Paclet Repository.
https://resources.wolframcloud.com/PacletRepository/resources/Wolfram/QuantumFramework/tutorial/QuantumOptimization.html

“Symbolic methods in quantum optimization: introduction to QAOA-in-QAOA”. Rodriguez, S.​ (2025). Wolfram Community, STAFF PICKS.
https://community.wolfram.com/groups/-/m/t/3462606

“Exploring the quantum approximate optimization algorithm (QAOA) through the max-cut problem​”. Rodriguez, S. (2025). Wolfram Community, STAFF PICKS.
https://community.wolfram.com/groups/-/m/t/3449705

Zhou, Z., Du, Y., Tian, X., & Tao, D. (2023). QAOA-in-QAOA: Solving Large-Scale MaxCut Problems on Small Quantum Machines. Physical Review Applied, 19(2).
https://doi.org/10.1103/physrevapplied.19.024027

Esposito, A., & Danzig, T. (2024). Hybrid Classical-Quantum Simulation of MaxCut using QAOA-in-QAOA. In 2024 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW) (pp. 1088–1094).
https://doi.org/10.1109/ipdpsw63119.2024.00180

Blekos, K., Brand, D., Ceschini, A., Chou, C.-H., Li, R.-H., Pandya, K., & Summer, A. (2024). “A review on Quantum Approximate Optimization Algorithm and its variants”. Physics Reports (Vol. 1068, pp. 1–66).
https://doi.org/10.1016/j.physrep.2024.03.002

Ruslan Shaydulin (2021). Tutorial on Combinatorial Optimization on Quantum Computers. YouTube.
https://www.youtube.com/watch?v=5bSH1JIqyko