numericaleft / feynmandiagram.jl

Computational graph representation of multiloop Feynman diagrams

FeynmanDiagram

FeynmanDiagram.jl is a Julia package designed to efficiently encode Feynman diagrams --- essential elements of Quantum Field Theory (QFT) —-- into compact computational graphs for fast evaluation. It employs Taylor-mode Automatic Differentiation (AD) specifically to implement field-theoretic renormalization schemes, a pivotal technique in QFT that significantly improves the convergence of Feynman diagrammatic series. This approach underscores the synergy between QFT and AI technologies, effectively addressing the sophisticated computational challenges in QFT.

Key Features

• Computational Graphs for Feynman Diagrams: Utilizes computational graphs analogous to those in Machine Learning (ML) architectures, where nodes represent mathematical operations and edges denote the tensor flow.

• Compiler Architecture: Implements a three-stage compiler process, analogous to modern programming language compilers, to transform Feynman diagrams into executable code, significantly enhancing the adaptability and computational efficiency of QFT calculations.

• Interdisciplinary Approach: Bridges QFT and AI tech stack by adapting ML algorithms, such as Taylor-mode AD, for QFT calculations, enhancing computational efficiency with tools like JAX, TensorFlow, PyTorch, and Mindspore.

Compiler Architecture Overview

In general, Feynman diagrams represent high-order integrals. The integrands are propagators and interactions composed with basic arithmetic operations (multiplication, addition, power, etc). The sequence of calculating the integrand by combining propagators and interactions with arithmetic operations can be represented as an algebraic computational graph. In this sense, the computational graph serves as an intermediate representation that standardizes the treatment of various diagram types, ensuring a consistent approach across different QFT calculations.

Drawing from these insights, the architecture of our compiler is intentionally crafted to process Feynman diagrams via a strategic, three-stage process that reflects the advanced design principles of the LLVM compiler architecture. This approach enhances the compiler's flexibility for a wide array of QFT computations and significantly improves its computational efficiency. The workflow is organized into three critical stages:

• Front End: In this initial phase, Feynman diagrams are generated and transformed into a standardized intermediate representation, taking the form of static computational graphs. This stage features two key algorithms for generating generic Feynman diagrams for weak-coupling expansions:

  • The Parquet module systematically organizes higher-order Feynman diagrams into a concise hierarchical structure of sub-diagrams, enabling efficient evaluation of repeated sub-diagrams. This module introduces algorithms for constructing streamlined computational graphs for two-, three-, and four-point vertex functions by exploiting the perturbative representations of the Dyson-Schwinger and parquet equations. This approach facilitates the analysis of a wide array of observables in quantum many-body problems.
  • The GV module, specifically aimed at many-electron problems with Coulomb interactions, incorporates the algorithm proposed in Nat Commun 10, 3725 (2019).
  • The front end also allows for the integration of new diagram types by users, enhancing its adaptability.

• Intermediate Representation: At this stage, the compiler applies optimizations and Automatic Differentiation (AD) to the static computational graph. This process is geared towards refining the graph for thorough analysis and optimization. The optimizations are aimed at streamlining the graph by removing redundant elements, flattening nested chains, and eliminating zero-valued nodes, among other strategies. The incorporation of Taylor-mode AD is critical for efficiently calculating high-order derivatives, essential for deriving diagrams for specific heat, RG flow equations, or renormalized Feynman diagrams.

• Back End: This final phase is responsible for translating the optimized graph into executable code that is compatible with a broad spectrum of computing environments. It supports various programming languages and facilitates seamless integration with different software and hardware ecosystems, significantly extending the compiler's utility across multiple platforms.

Usage

Installation

To install the package, you can simply use the following command in the Julia package manager:

using Pkg
Pkg.add("FeynmanDiagram")

Example: Self-energy in FrontEnds

The algorithms in FrontEnds generate Feynman diagrams for weak coupling expansions, supporting various diagram types, such as self-energy, polarization, 3-point vertex function, and 4-point vertex function diagrams. The internal degrees of freedom can be either the loop variables (e.g., momentum or frequency) or the site variables (e.g., imaginary-time or lattice site).

The example below demonstrates how to generate all two-loop self-energy diagrams by Parquet and optimize their computational graphs.

Diagram generation with the Parquet algorithm

using FeynmanDiagram
import FeynmanDiagram.FrontEnds: NoHartree

# Define a parameter for two-loop self-energy diagrams in the momentum and the imaginary-time representation. Exclude any diagrams containing Hartree subdiagrams. 
para = Parquet.DiagPara(type = Parquet.SigmaDiag, innerLoopNum = 2, hasTau = true, filter=[NoHartree,]);

# Construct Feynman diagrams within a DataFrame utilizing the parquet algorithm. The resulting sigmadf DataFrame comprises two components: the instantaneous part and the dynamic part of the self-energy.
sigmadf = Parquet.build(para) 
2×4 DataFrame
 Row │ diagram                            extT    hash   type
     │ Graph…                             Tuple…  Int64  Analytic…
─────┼─────────────────────────────────────────────────────────────
   1 │ 18721-ΣIns, k[1.0, 0.0, 0.0], t(…  (1, 1)  18721  Instant
   2 │ 18722-ΣDyn, k[1.0, 0.0, 0.0], t(…  (1, 2)  18722  Dynamic

# Optimize the Graph for the given Feynman diagrams.
optimize!(sigmadf.diagram); 

Construct renormalized Feynman diagrams using Taylor-mode AD

The example code below demonstrates how to build renormalized Feynman diagrams for the self-energy including Green's function and interaction counterterms using Taylor-mode AD.

# Set the renormalization orders. The first element is the maximum order of the Green's function counterterms, and the second element is the maximum order of the interaction counterterms.
renormalization_orders = [2, 1];

# Define functions that determine how the differentiation variables depend on the properties of the leaves in your graphs, identifying `BareGreenId` and `BareInteractionId` properties as the Green's function and interaction counterterms, respectively.
leaf_dep_funcs = [pr -> pr isa FrontEnds.BareGreenId, pr -> pr isa FrontEnds.BareInteractionId];

# Generate the Dict of Graph for the renormalized self-energy diagrams with the Green's function counterterms and the interaction counterterms.
dict_sigma = taylorAD(sigmadf.diagram, renormalization_orders, leaf_dep_funcs)
Dict{Vector{Int64}, Vector{Graph}} with 6 entries:
  [0, 0] => [18822Ins, k[1.0, 0.0, 0.0], t(1, 1)=0.0=Ⓧ , 19019Dyn, k[1.0, 0.0, 0.0], t(1, 2)=0.0=⨁ ]
  [2, 1] => [18823Ins, k[1.0, 0.0, 0.0], t(1, 1)[2, 1]=0.0=Ⓧ , 19020Dyn, k[1.0, 0.0, 0.0], t(1, 2)[2, 1]=0.0=⨁ ]
  [2, 0] => [18824Ins, k[1.0, 0.0, 0.0], t(1, 1)[2]=0.0=Ⓧ , 19021Dyn, k[1.0, 0.0, 0.0], t(1, 2)[2]=0.0=⨁ ]
  [1, 1] => [18825Ins, k[1.0, 0.0, 0.0], t(1, 1)[1, 1]=0.0=Ⓧ , 19022Dyn, k[1.0, 0.0, 0.0], t(1, 2)[1, 1]=0.0=⨁ ]
  [1, 0] => [18826Ins, k[1.0, 0.0, 0.0], t(1, 1)[1]=0.0=Ⓧ , 19023Dyn, k[1.0, 0.0, 0.0], t(1, 2)[1]=0.0=⨁ ]
  [0, 1] => [18827Ins, k[1.0, 0.0, 0.0], t(1, 1)[0, 1]=0.0=Ⓧ , 19024Dyn, k[1.0, 0.0, 0.0], t(1, 2)[0, 1]=0.0=⨁ ]

Example: Compile Feynman diagrams to different programming languages

The Back End architecture enables the compiler to output source code in a range of other programming languages and machine learning frameworks. The example code below demonstrates how to use the Compilers to generate the source code for the self-energy diagrams in Julia, C, and Python.

# Access the self-energy data for the configuration with 2nd-order Green's function counterterms and 1st-order interaction counterterms.
g_o21 = dict_sigma[[2,1]]; 

# Compile the selected self-energy into a Julia RuntimeGeneratedFunction `func` and a `leafmap`.
# The `leafmap` associates vector indices of leaf values with their corresponding leaves (propagators and interactions). 
func, leafmap = Compilers.compile(g_o21);

# Export the self-energy's source code to a Julia file.
Compilers.compile_Julia(g_o21, "func_o21.jl");

# Export the self-energy's source code to a C file.
Compilers.compile_C(g_o21, "func_o21.c");

# Export the self-energy's source code to a Python file.
Compilers.compile_Python(g_o21, "func_o21.py");

Computational Graph visualization

Tree-type visualization using ETE3

To visualize tree-type structures of the self-energy in the Parquet example, install the ETE3 Python package, a powerful toolkit for tree visualizations.

Execute the following command in Julia for tree-type visualization of the self-energy generated in the above Parquet example:

plot_tree(sigmadf.diagram, depth = 3)

For installation instructions on using ETE3 with PyCall.jl, please refer to the PyCall.jl documentation on how to configure and use external Python packages within Julia.

Graph visualization using DOT

The Back-End's Compilers module includes functionality to translate computational graphs into .dot files, which can be visualized using the dot command from the Graphviz software package. Below is an example code snippet that illustrates how to visualize the computational graph of the self-energy from the previously mentioned Parquet example.

# Convert the self-energy graph from the `Parquet` example into a dot file.
Compilers.compile_dot(sigmadf.diagram, "sigma_o2.dot");

# Use Graphviz's dot command to create an SVG visualization of the computational graph.
run(`dot -Tsvg sigma_o2.dot -o sigma_o2.svg`);

The resulting computational graphs are depicted as directed acyclic graphs by the dot compiler. In these visualizations, the leaves are indicated by green oval nodes, while the intermediate nodes take on a rectangular shape. On the penultimate level of the graph, the left blue node with the Sum operator signifies the dynamic part of the self-energy, and the right blue node denotes the instantaneous part of the self-energy.

License

FeynmanDiagram.jl is open-source, available under the MIT License. For more details, see the license.md file in the repository.

feynmandiagram.jl's Issues

Optimization for Computational graph

After I visualized the computational graph through the compile_dot function, by inspecting the simplest graph (1,0,0) for GV and parquet algorithm, I found some optimizations don't work. As the graph show below:
For the parquet graph:

It's obvious to see that there are still trivial unary chains,
The same situation happen for the GV graph

The reason may be that the factor of some nodes is not 1.

Add essential IR tree operations

Write unit tests and functions for essential tree transformations in the intermediate (computational graph) representation:

Transformations

Transformations for tree sanitization:

• prune_unary — Prune trivial unary operations, e.g., $(+a) \rightarrow a$ and $(*a) \equiv (1*a) \rightarrow a$.
• inplace_prod — Propagate subgraph factors up multiplicative chains and merge.
• merge_prefactors — Factorize multiplicative graph prefactors, e.g., $a * G + b * G \rightarrow (a + b) * G$.
• merge_chains — Merge operator chains into a single operation, e.g., $a+(+(\cdots+(+b)\cdots)) \rightarrow a+b$ and $G_1 +(+(\cdots+)+)\; G_2 \rightarrow G_1 + G_2$.
• merge_operators — Use the associativity of an operator $\mathcal{O}$ to "deparenthesize" an expression, i.e., join two sub-operations into a single $n$-ary form. Use with $\mathcal{O} = \bigotimes$ requires Base.:*(g1, g2).
• split_operator — Use the associativity of an operator $\mathcal{O}$ to "parenthesize" an expression, i.e., split the $n$-ary form of a node into two sub-operations. Use with $\mathcal{O} = \bigotimes$ requires Base.:/(g1, g2).

---

High-priority transformations for the Parquet front-end:

• relabel — Update the label(s) of quantum operators in a graph and its subgraphs, e.g., {{1,2}=>3, 3=>4}.
• standarize_labels — Finds all labels in a graph and its subgraphs (e.g., 1,4,5,7...) and then relabels them in standard order 1,2,3,4,...
• unghost_legs — A special case of split_operator for graph decomposition $(\mathcal{O} = \bigotimes)$ which does not require Base.:/(g1, g2). Extracts a vertex with ghost legs and its connected propagators into a new $\bigotimes$ node, thereby "unghosting" its legs.
• replace_subgraph — Replace one subgraph of a graph g with another one of possibly different topology.

---

Low-priority transformations requiring associative graph (de)composition to be defined via Base.:*(g1, g2) and Base.:/(g1, g2):

• expand — Expand a factorized graph expression: $(G_1 + G_2) * W \mapsto G_1 * W + G_2 * W$. The recursive depth should be specified by an argument depth which is 1 by default.
• factorize — Factorize a graph expression (the inverse of expand): $G_1 * W + G_2 * W \mapsto (G_1 + G_2) * W$. The recursive depth should be specified by an argument depth which is 1 by default. Requires Base.:/(g1, g2).
• tofeynmanrep — Convert a general graph to the Feynman representation by chaining expand and merge_chains.
• fromfeynmanrep — Optimize a graph in the Feynman representation by chaining factorize and merge_chains.

Definitions/Examples

---

---

Add quantum operators

Implement a struct QuantumOperator for creation/annihilation operators. Further use this struct to define composite operators and operator pairings via an algebra for QuantumOperator instances.

The aims are (1) to rigorously define generic internal/external vertices and (2) to support derivation of ExternalVertex types for a diagram given user-supplied Feynman rules (a list of valid internal vertices for the problem).

rename loop to K, and Tau to X

Some suggestions:

1. Loop --> K
2. Tau ---> X

Feature wishlist for new API

Frontend

A user-extensible interface to the IR.

• Inline parser (deserializer)
• Serializer (tostring method)
• Batch parser for TOML config file(s)
• Generic support for n-particle reducibility

Intermediate Representation (IR)

A computational tree representation for generic Feynman diagram(s).

Tree operations

• Auto-differentiation (e.g., RG flows)
• Auto-integration (e.g., Matsubara summation for Feynman diagrams and/or general IR trees)

Tree transformations

• Leaf replacement (subdiagram expansion)
• Flattening transformation: reduce tree to n layers
• (Un)distribute transformation: move addition (multiplication) to subdiagram root
• Filtering transformation: diagram removal based on reducibility rules / custom filter
• IR to Feynman diagram representation (special case of flatten + distribute transformation)

Backend

Optimization of IR trees and compilation to expression trees. For optimizations, we can either follow or interface with aesara.

Tree optimizations

• Remove duplicate nodes
• Merge (combine redundant operations)
• GEMM optimization for matrix addition/multiplication
• Support for matrix product states for higher-order vertex functions

GraphVector struct

Implement a data structure to handle vector graphs:

NOTE: All graphs in the vector should have the same number of external legs, but the labels of the external legs can be different.

• Add the following struct
  struct GraphVector{F, W} <: AbstractVector
  graphs::Vector{Graphs{F, W}}
  end
• Implement vector indexing through julia iterator API
• A function to group the graphs with the same external legs with a given index. For example,

For a GraphVector, gv = GraphVector([g(1, 2), g(1, 3), g(2, 3)]). We expect the following:

1. merge(gv, indices = [1, ]) = [GraphVector([g(1, 2), g(1, 3)]), GraphVector([g(2,3), ])]
2. merge(gv, indices = [2, ]) = [GraphVector([g(1, 2), ]), GraphVector([g(2,3), g(1, 3)])]
3. merge(gv, indices = [1, 2]) = [GraphVector([g(1, 2), ]), GraphVector([g(2,3), ]), GraphVector([g(1, 3), ])]

• construct Feynman diagram from a set of GraphVectors.

function feynman_diagram(vertices::Vector{GraphVector}, topology::Vector{Vector{Int}};
external=[], factor=one(_dtype.factor), weight=zero(_dtype.weight), name="", type=:generic)

Compile hand-drawing Feynman diagram into diagram tree

In the src/diagramBuilder/fromFile, add code to compile the user-defined Feynman diagram into a diagram tree.

The diagram file could be generated from the external package or manually created by the user. The code reads the file and compiles it into a diagram tree.

Question on the function `merge_prefactor` and `Base.+`

I found that the function merge_prefactor (in #58) has to confronted the different cases below and it seems a little messy, for g=propagator(𝑓⁺(1)𝑓⁻(2), factor=1):
(a)k1=(1*g+2*g)
(b)k2= linear_combination(g,g,1,2)
The structures of them are different:
(1) For k1:

(2) For k2:

I suggest these two cases should be unified to the second one in the Base.+ function.

Update transformations to support generic computational graphs

Many transformations in transform.jl and related optimizations in optimize.jl apply only to expression trees, i.e., they assume each node has only one parent. Functions relabel / relabel! and standardize_labels / standardize_labels! are already sufficiently general. All other functions must be corrected to support generic computational graphs, which may include nodes with multiple parents.

feynman_diagram constructor from vector of graphs

implement a function to generate a feynman_diagram from a vector of graphs with the following signature.

function feynman_diagram(vertices::Vector{Graph}, topology::Vector{Vector{Int}};
external=[], factor=one(_dtype.factor), weight=zero(_dtype.weight), name="", type=:generic)

vertex3 from the parquet equation may need an additional sign

The one-loop vertex3 calculated from the following source code gives F = -0.5123(18). Since z factor is expected to be 0.662, so that we expect F = 1/z-1=0.51057. It indicates an additional sign is needed for the fermionic vertex correction.

https://github.com/numericalEFT/FeynmanDiagram.jl/blob/81a7fa415eb90454057662f53d9d1ad613cf7cdb/example/vertex4/ver3_l.jl

add support for autodiff operation on the diagram tree

builder for self-energy, polarization and vertex3

Feynman rule convention

in the Feynman rule docs, https://numericaleft.github.io/FeynmanDiagram.jl/dev/manual/feynman_rule/, we are currently using G^{-1}=g^{-1}+\Sigma. However, in many books, it should be G^{-1}=g^{-1}-\Sigma in books.

More generic diagram tree

1. In general, the variables of the diagram object could be the combination of the following:

• Momentum
• Frequency
• Pair of space
• Pair of time
  The variables could be internal or external.

2. Generic diagram components:

• Propagator
• Vertex with N legs (note that one vertex may have multiple values, e.g., W has an instantaneous interaction and a retarded one)
• Nodes (a sum/product of propagators and vertices).

Note that when evaluating a node, one may need to evaluate an additional weight factor based on the property of a node.

One may implement the above codes in the following roadmap:

• An abstract struct called Pool to host an array of unique objects (could be a variable, a propagator, or a vertex). Implement the pool management operations (add an object, remove an object, find an object, ...)
• Pool struct of the propagators.
• Pool struct of the vertices.
• Diagram tree based on the propagator and vertex pools.
• Diagram tree evaluation.

More efficient algorithm to evaluate diagram tree

In the naive algorithm, the entire diagram tree is recalculated. In a more advanced algorithm, one should only evaluate part of the tree in each MC update.

Basic IR transform required by the Parquet Algorithm

[ ] relabel(g::Graph, map::DIct{Int, Int}):
this function maps the labels of the quantum operators in g and its subgraphs to other labels. For example, map = {1=>2, 3=>2} will find all quantum operators with labels 1 and 3, and then map them to 2.

[ ] standarize_labels(g::Graph):
this function first finds all labels involved in g and its subgraphs (for example, 1, 4, 5, 7, ...), then relabel them in the order 1, 2, 3, 4, ....

Parquet algorithm to Diagram Tree

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Error when generating (3,0,0) vertex4 diagram

using ElectronLiquid, FeynmanDiagram
para = UEG.ParaMC(rs=0.5, beta=25.0, Fs=-0.0, order=4, mass2=3.0, isDynamic=false, dim=2);
diagram = Ver4.diagram(para, [(3,0,0)]; channel=[PHr, PHEr, PPr], filter=[NoHartree, Proper])

ERROR: DimensionMismatch: dimensions must match: a has dims (Base.OneTo(5),), b has dims (Base.OneTo(16),), mismatch at 1

refactor DiagramID

The current implementation of DiagramID duplicates with DIagPara.type.

DiagramID should have larger scope than DiagPara.type. To be more specific,

DiagPara.type --> Ver4, Ver3, Green, Interaction, Sigma, ...
DiagramID --> N-body Vertex/Green/connected-Green in momentum-time representation, spacetime-representation, ...

So it is probably better to call it Representation. An example is the following,

#generic Veretx function in KX representation
struct VerKX
    para::DiagPara #DiagPara don't need weightType
    spin::Vector{Int} #[Up, Up], [Up, Down], ...
    channel::Channel # particle-hole, particle-hole exchange, particle-particle, irreducible, and their counterterms
    extK::Vector{Vector{Float64}}
    extT::Vector{Int}
    order::Vector{Int}   #[2, 0, 1] means loop-order 2, and one derivative of the second type.
end

it should come with some help functions (for example, a function to check if the vertex is instant or not).

The Diagram struct should be defined as,

struct Diagram{W, ID<:DiagramId}
      hash::Int
      name::Symbol
      id::ID
      operator::Operator
      factor::W
      subdiagram::Vector{Diagram{W, ID}}

      weight::W

Need some help functions (for example, isbare(diag) = (length(diag.subdiagram)==0))

add support for proper diagram and one-particle-reducible diagram