# Combining deep learning and street view imagery to map smallholder crop
types

Jordi Laguarta1, Thomas Friedel2, Sherrie Wang1

###### Abstract

Accurate crop type maps are an essential source of information for monitoring
yield progress at scale, projecting global crop production, and planning
effective policies. To date, however, crop type maps remain challenging to
create in low- and middle-income countries due to a lack of ground truth
labels for training machine learning models. Field surveys are the gold
standard in terms of accuracy but require an often-prohibitively large amount
of time, money, and statistical capacity. In recent years, street-level
imagery, such as Google Street View, KartaView, and Mapillary, has become
available around the world. Such imagery contains rich information about crop
types grown at particular locations and times. In this work, we develop an
automated system to generate crop type ground references using deep learning
and Google Street View imagery. The method efficiently curates a set of street
view images containing crop fields, trains a model to predict crop type by
utilizing weakly-labelled images from disparate out-of-domain sources, and
combines predicted labels with remote sensing time series to create a wall-to-
wall crop type map. We show that, in Thailand, the resulting country-wide map
of rice, cassava, maize, and sugarcane achieves an accuracy of 93%. As the
availability of roadside imagery expands, our pipeline provides a way to map
crop types at scale around the globe, especially in underserved smallholder
regions.

##  1 Introduction & Background

Figure 1: Top: Street-view images of roadside occlusions present between the
car-mounted camera and fields. Bottom: Street-view images after the automated
filtering process for the four major crop types in Thailand. Figure 2:
Overview of the methods presented in this paper to create a Thailand-wide crop
type map. Example field points, ground reference labels, and crop type map are
shown for the district of Sawang Ha.

Ensuring global food security is one of the major challenges we will face this
century, especially in the face of a changing climate and a growing climate
population (Becker-Reshef et al. 2023; Rezaei et al. 2021). Accurate crop type
maps are an essential source of information for monitoring yield progress
(Becker-Reshef et al. 2023), projecting global crop production (Lobell and
Gourdji 2012), and planning effective policies (Haasnoot et al. 2013).
However, only a handful of countries—mostly high-income—have had the budget to
collect large-scale ground data and develop crop type maps (Han et al. 2012;
Fisette et al. 2013; Cantelaube and Carles 2014). Meanwhile, regions with
smallholder farms, which provide a living for two-thirds of the world’s rural
population of 3 billion (Lowder, Skoet, and Raney 2016) and produce 80% of the
world’s food (Economic et al. 2014), continue to lack such maps. The majority
of smallholder farms are located in middle- and low-income countries, where
expensive ground data on crop types remains scarce (Tseng et al. 2021; Wang et
al. 2020; Lee et al. 2022).

To address the high cost of acquiring ground reference labels, the remote
sensing and machine learning communities have started to explore non-
traditional sources of crop type data. One such source is roadside images
through services like Google Street View (GSV), Bing Maps StreetSide,
Mapillary, KartaView, Tencent Street View, and Baidu Total View. The images
are captured by dash cams or panoramic cameras mounted on cars; depending on
the service, they are crowdsourced or collected by dedicated fleets. Today,
they are low-cost to access, available in almost every country, and updated
every few years.

Recent works have used roadside imagery in applications ranging from urban
morphology to real estate to air quality prediction (Biljecki and Ito 2021).
Most relevant to our work, Paliyam et al. (2021) deployed cameras mounted on
the hoods of vehicles in Africa, Wu et al. (2021) used smartphones mounted on
car windows, and Yan and Ryu (2021) used GSV to create crop type ground
reference labels. Recently, the European Space Agency and WorldCereal
Consortium released a crop type map for cereals leveraging GSV to manually
create a validation set (Van Tricht et al. 2023).

Figure 3: Spatial and temporal distribution of Google Street View in Thailand.
Left: Hexbin plot of GSV availability across Thailand. The zoomed-in panel
shows the location of street-view images overlaid on a satellite basemap in
the district of Sawaeng Ha. Right: Availability of street-view images in
Thailand by month, with a clear rise in availability since 2022 and a total of
over 3 million images. During the wet season (May–October) shown in the blue
box, 1.5 million images are available.

However, existing approaches remain either small in scale (i.e., rely on
manual labeling) or are difficult to deploy in smallholder regions. Challenges
in smallholder regions include: complex road networks in rural areas compared
to the grid system present in the US Midwest (Yan and Ryu 2021) and vegetation
and man-made occlusions blocking the view between the road and fields (Figure
1). Furthermore, as different regions grow different crops, current methods
require hand-labelling a new crop type roadside dataset when encountering a
new region.

In this paper, we propose a cost-effective automated deep learning pipeline to
generate crop type references and remote sensing crop type maps with minimal
manual labeling. We create a method to auto-generate field coordinates at
scale, and then filter GSV images to create a large dataset of geo-tagged
field images (Section 2.2). Next, a scraper automatically curates a training
set of weakly-labeled images from disparate out-of-domain sources to train a
CNN crop type discriminator on street-view images for the crop types of
interest (Section 2.3). The crop type of field images is inferred by the CNN
to generate ground referenced labels as a training set for a remote sensing
crop type mapper (Sections 2.3-2.4). Once trained, the remote sensing model
outputs crop type maps (Section 2.4).

This work shows for the first time smallholder crop type mapping on a country
scale using street-level imagery. We tested our approach in Thailand for the
May 2022 to October 2022 wet growing season on the region’s four major crops
and created a ground truth set of 1600 GSV images labelled by a plant taxonomy
expert. A total of 81,000 ground reference points were generated using our
deep learning pipeline to train the remote sensing crop type model on the
whole country. The crop type maps achieved an overall accuracy of 0.93 and an
F1 score of 0.92. The approach is orders of magnitude cheaper than traditional
survey-based methods, scalable, high accuracy, and automated — with expert
hand-labeling only necessary to create a test set. The approach extends to
non-GSV street-level imagery as well.

We open source all code, datasets, and crop type maps generated here: [link
omitted for blind review].

##  2 Datasets & Methods

In brief, our method generates proposal locations along streets, curates a set
of images containing crop fields (Section 2.2), predicts the crop type within
the street-level image with minimal manual labeling (Section 2.3), and
combines these labels with remote sensing time series to create a wall-to-wall
crop type map (Section 2.4) (Figure 2).

###  2.1 Study area

We chose Thailand as the area of study, because it is simultaneously dominated
by small-scale farms and has a high availability of Google Street View images.
GSV images in Thailand grew 700% from March 2022 to May 2023 to a total of
3.9M images, of which 2.8M are in cropland areas across the country (Figure
3). Two-thirds of cropland in the country is used to grow rice, while
sugarcane, cassava, and maize are the next most abundant crops at just under
10% of crop area each (Figure 1) (FAO 2021). Rice is grown in two seasons, wet
and dry; we limit this work to the wet season.

###  2.2 Finding street-level images of crop fields

#### Extrapolate equidistant street points

The first step in the pipeline is to gather the latitude and longitude
coordinates of candidate fields and their corresponding roadside image along
all roads in the country. We started with Open Street Map (OSM), a worldwide
open geographic database with high coverage and completeness (OpenStreetMap
contributors 2017), and used Overpass API to query all OSM ways within
Thailand. In order to maximize recall of GSV field images, we generated
equidistant points at 10m steps along OSM ways. Examples of points sampled
from OSM street nodes in a chosen area are shown in Figure A.1.

#### Filter points using land cover map

Next, we used existing land cover maps to filter for street points near
farmland and remove field points that are not visible from the street due to
obstructions (Figure 1). In particular, we used the European Space Agency’s
WorldCover 10m v100, which classifies global land cover at 10 meter resolution
for the year 2020. The map contains 11 land cover classes, including tree
cover, shrubland, grassland, built-up, and cropland. We accessed WorldCover
and filtered candidate points using functions available in Google Earth
Engine.

Starting from the OSM-derived equidistant points, we removed points containing
no cropland within a 10m radius. An inspection of 200 remaining points
revealed that, in 30% of GSV images, crop fields were still blocked by trees
and other vegetation next to the road. We therefore also removed points
containing any tree cover within the same 10m radius.

Finally, we are only interested in roadside images taken during the growing
season, when crop types may be visible from the road. For this work, only GSV
images captured during the wet growing season in Thailand (May 2022–October
2022) were considered. The above approach can be modified to any region and
growing season worldwide by changing the date range and shapefile of the
region of interest.

#### Extrapolate camera heading and field point

Due to the grid layout of streets in the US Midwest, prior work using GSV for
crop type mapping obtained roadside images facing crop fields and points
within crop fields by extrapolating street points in north/south and east/west
directions (Yan and Ryu 2021). Since street layouts are more complex in
smallholder regions, we first found the street bearing $\theta$ using the
haversine formula (see Section A.2), and then computed the heading for the
camera to face the two fields on either side of the street as $\theta\pm 90$
degrees (Figure A.1). We empirically determined that 30m was the best distance
to extrapolate street points to crop field interior points (Table A.2).

Figure 4: Schematic of the process to generate ground reference labels. Each ground reference is composed of crop type and geocoordinates from street-view images. | Percent of Dataset | Number of Samples
---|---|---
Crop Type | Planted Area | Expert Street View | WebCC | iNaturalist | Expert Street View | WebCC | iNaturalist
Rice | 67% | 79% | 24% | 16% | 1261 | 659 | 1396
Sugarcane | 7% | 9% | 25% | 22% | 144 | 679 | 1882
Cassava | 8% | 5% | 21% | 42% | 81 | 584 | 3662
Maize | 8% | 7% | 30% | 20% | 111 | 832 | 1698
Total | 90% | 100% | 100% | 100% | 1597 | 2754 | 8638
Table 1: Crop type distribution for the various datasets used to train and
test a classifier on street-view images in Thailand. Planted area is obtained
from annual national statistics released by the FAO. Street View refers to GSV
images randomly sampled from Thailand and manually labelled with crop type.
WebCC are images scraped from the web. iNaturalist are images selected from an
online biodiversity database.

#### Classify in-field images

After filtering out points near trees and finding the appropriate camera
heading, we found that 58% of GSV images still had an obstructed view of a
field due to small bushes between the road and field not detected by satellite
land cover maps (Table 3). We therefore labelled 2986 images into {field, not
field} and trained a binary classifier using a ResNet-18 pre-trained on
ImageNet. Hyperparameters used include an Adam optimizer and learning rate of
$0.001$ for 15 epochs. The model classified field images from non-field images
with 95% recall and 98% precision. Candidate street-view images were
downloaded ($n=224,000$) and classified by the field/not-field CNN. Those
labelled as field ($n=89,000$) were used as input in Section 2.3 to be
classified into crop types.

###  2.3 Predicting crop type in street-level images

##### Compile training set from the web

Data annotation to train a new classifier requires manual labeling, which is
time consuming and, for crop type classification, requires domain expertise.
Fortunately, the internet can serve as a source to rapidly obtain a training
set of images that is large, low-cost, and diverse with real-world settings
(Krause et al. 2016).

We compiled training sets of crop type images from two sources: Google Images
Creative Commons (hereafter “WebCC”) and iNaturalist. For each crop type,
images were queried in Google Images by searching for the crop name followed
by “field” (e.g., “rice field”). Returned images were labeled with the queried
crop type. Meanwhile, iNaturalist is an online community of naturalists and
citizen scientists who contribute photos of biodiversity across the globe. Its
database contains over 161 million observations. While not targeted toward
agriculture, iNaturalist contains images of crops, which can be downloaded via
an API (Huerta-Ramos and Luštrik 2021). In total, thousands of image-label
pairs were collected from the two platforms (Table 1).

Although we tried to compile training images that look similar to GSV, online
images varied in quality, size, and relevance. We observed many images of
people with fields in the background, people holding crops, and label noise in
WebCC (Figure 5). Fortunately, CNNs tend to be robust to noise and can perform
well if enough signal from the desired task is present during training (Krause
et al. 2016; Carlini and Wagner 2017; Szegedy et al. 2013; Wei et al. 2021).
Therefore, we did not manually clean the dataset.

Figure 5: Example images from the two online datasets. Although some images
are of crop fields similar to the target street-view task, many images are
either close-ups of a plant (especially in iNaturalist), a single plant
instead of a field, or label noise (especially in WebCC).

##### Manually label crop type ground references

We collected a test set of GSV images randomly sampled by geography and had
them manually labeled by a plant taxonomy expert. To avoid data leakage in the
remote sensing evaluation from two points belonging to the same field, we
ensured all field points were more than 100m away from each other. We labelled
1922 images into the following classes: cassava, maize, rice, sugarcane,
unknown, non-crop, unsupported crop, and additional. The expert-labeled
dataset served two purposes, (1) to compare performance of the CNN street-view
classifier across different training sets and (2) to evaluate the remote
sensing crop type classifier.

##### Train a street-level crop type classifier

We trained models on 5 different datasets, 3 of which contained only online
images automatically labeled by the corresponding search label. We use WebCC
and iNaturalist to denote models trained on Google Images results and
iNaturalist images, respectively. iNaturalist + WebCC was trained by combining
the WebCC and iNaturalist datasets. Expert labeled was trained on 60% of the
expert-labeled GSV test set. Finally, Combined merged the same 60% set of GSV
with WebCC and iNaturalist. For all models, the remaining 40% of the expert-
labeled GSV dataset was split evenly into validation and test sets.

We used data augmentation techniques to accommodate images that varied in size
across datasets and help the model focus on crop type features rather than a
crop’s position, zoom, or orientation in an image. During training, all images
were expanded to $600\times 600$ px and horizontally flipped at random,
followed by a random crop of $300\times 300$ px.

We selected a ResNet-50 architecture pre-trained on ImageNet for ease of
reproducibility. All models were trained on 4-class classification with cross-
entropy loss for 15 epochs, with a learning rate starting at $0.001$. A cosine
annealing learning rate scheduler was used to dynamically adjust the learning
rate and help improve model convergence.

#### Predict crop type in street-view images

Trained models were used to classify street-view images ($n=89,000$) into crop
types. We took sliding windows across each image, with a window of $300\times
300$ pixels and a stride of $50$. We used the mode of high probabilities (MHP)
to select the predicted class from all the sliding windows in an image. For
each window, if p is the vector of softmax probabilities and $c$ is the crop
type with the highest probability, then the window casts one vote for class
$c$ if $\textbf{p}_{c}$ exceeds some threshold $\tau$. If none of its
probabilities exceed $\tau$, then the window casts no votes. The final
prediction for the image is the crop type with the most votes across all
windows.

Each clean street-view image was classified by the 5 models with sliding
window and MHP and paired with its field coordinate to generate the ground
reference training sets.

Pipeline Step (Section) | Dataset Size
---|---
Candidate field points (2.2) | 98 million
Filtered field points (2.2) | 3.9 million
GSV field images available (2.2) | 2.8 million
Street-view images downloaded (2.2) | 224,000
Clean street-view images (2.2) | 89,000
Labelled field points (2.3) | 89,000
MHP thresholding (2.3) | 81,000
Table 2: Summary of dataset size through the point and image filtering steps
in the pipeline. Starting with a road network of Thailand and ending with
81,000 ground reference labels for the four major crop types.

###  2.4 Remote sensing-based crop type mapping

#### Feature extraction via harmonic regression

Consistent with existing state-of-the-art methods (Inglada et al. 2015; Jin et
al. 2019), we used Sentinel-2 satellite time series for country-wide crop type
classification. Sentinel-2 images capture 13 optical bands at 10-60m
resolution on a 5-day cycle, and their time series contain information on crop
phenology that allows different crop types to be distinguished.

Method | Precision | Recall
---|---|---
Baseline (Yan and Ryu 2021) | 0.07 | 0.25
Camera heading | 0.14 | 1.00
Cropland filter | 0.31 | 1.00
Cropland and tree cover filters | 0.42 | 0.99
All + field classifier | 0.98 | 0.95
Table 3: Precision and recall for filtering out non-field GSV images. Baseline is our implementation of the approach in Yan and Ryu (2021). Each Method includes the methods above. All methods previous to the field classifier are filters applied before downloading the GSV image. | Street-view Test Set Metrics
---|---
Training Dataset | Overall Acc | Overall F1 | Cassava F1 | Maize F1 | Rice F1 | Sugarcane F1
Baseline: Most common class | 0.703 | 0.696 | 0.000 | 0.000 | 0.882 | 0.000
WebCC | 0.829 | 0.822 | 0.098 | 0.403 | 0.980 | 0.536
iNaturalist | 0.832 | 0.814 | 0.110 | 0.413 | 0.982 | 0.542
iNaturalist + WebCC | 0.938 | 0.960 | 0.885 | 0.707 | 0.992 | 0.797
Expert labeled | 0.939 | 0.963 | 0.813 | 0.882 | 0.978 | 0.854
Combined | 0.959 | 0.973 | 0.989 | 0.692 | 0.990 | 0.928
Table 4: Performance on crop type classification from street-view images for
the four major crops in Thailand. Models were trained on combinations of three
different training datasets: WebCC, iNaturalist, and Expert labeled GSV
images. The MHP threshold was set to 0.9.

We used Google Earth Engine to export Sentinel-2 L2A time series in Thailand
from May 1 – October 31, 2022. Four spectral bands were used: Red Edge 4, SWIR
1, SWIR 2, and NIR. We added the green chlorophyll vegetation index
($\text{GCVI}=\text{NIR}/\text{GREEN}-1$), as prior work showed it to be a
valuable feature for crop type classiﬁcation. We used the Cloud Probability
band to remove cloudy days before extracting time-series for each ground
reference (Figure A.3).

We used the harmonic regression to extract frequency-domain features from time
series of varying lengths for input to machine learning. Equivalent to the
discrete Fourier Transform, harmonic regression generates robust features for
crop type classification (Shumway, Stoffer, and Stoffer 2000; Jin et al. 2019;
Ghazaryan et al. 2018; Jakubauskas et al. 2001; Wang, Azzari, and Lobell 2019;
Wang et al. 2020). We applied a 3rd order harmonic regression to each band
independently to arrive at a total of 35 features.

#### Country-wide crop type classification

We trained random forests (with $500$ trees) on crop type classification, with
the harmonic coefficients as input and the crop type labels generated by the
street-view CNN as output. Random forests (Breiman 2001) are frequently used
in remote sensing applications for their high accuracy and computational
efficiency (Gislason, Benediktsson, and Sveinsson 2006; Azzari and Lobell
2017; Ok, Akar, and Gungor 2012).

The crop type ground references ($n=81,000$, Section 2.3) were used to train
the remote sensing crop type classifier. We trained 5 random forests, one for
each ground reference set (Table 4). A sixth model was trained only on the
expert-labeled dataset ($n=984$). Compared to automatically-labeled points,
expert-labeled points are more accurate but almost 100 times smaller in
quantity.

##  3 Results

###  3.1 Automatically generated field points and street-view images

We obtained 98 million points across Thailand by computing 10m equidistant
points from OSM ways. Of these points, 3.9 million occurred near cropland
without tree cover, and 2.8 million had GSV images available (Table 2). Of
these street-view images, we downloaded 224,000 and classified 89,000 into
images of fields and 135,000 into not fields. The 89,000 field-view images
were classified by crop type (rice, maize, cassava, sugarcane) and passed
though the MHP threshold to create 81,000 point ground reference labels. The
remaining 8,000 images were removed at MHP thresholding because none of their
sliding window predictions had a softmax probability greater than 0.9
indicating label uncertainty.

To understand how important land cover filtering, camera heading, and
field/not field classification were for finding GSV images of crops, we
calculated the precision and recall at various steps along the pipeline (Table
3). We also implemented the method developed by Yan and Ryu (2021) in the US
Midwest, which did not include these filtering steps, as a baseline. In a
sample of random GSV images in Thailand, we found that the baseline method
achieves a precision of 0.07 and recall of 0.25; in other words, 93% of
downloaded GSV images did not show a crop field, and 75% of images of crop
fields were missed. This illustrates the complexity of smallholder landscapes
compared to industrial agriculture. By contrast, using the correct camera
heading on points generated 10m apart from OSM ways yielded a precision of
0.14 and recall of 1.00. Filtering by cropland alone improved precision to
0.31, and removing points near trees further improved precision to 0.42.
Finally, training a CNN to classify field/not field improved precision to 0.98
and lowered recall slightly to 0.95 — a worthwhile trade-off.

###  3.2 Weakly supervised creation of crop type ground references

| Remote Sensing Test Set Metrics
---|---
Ground Ref. Train Dataset | Overall Acc | Overall F1 | Cassava F1 | Maize F1 | Rice F1 | Sugarcane F1
Baseline* | 0.695 | 0.456 | 0.000 | 0.000 | 0.753 | 0.000
Auto w/ WebCC† | 0.794 | 0.801 | 0.126 | 0.392 | 0.970 | 0.510
Auto w/ iNaturalist† | 0.783 | 0.796 | 0.132 | 0.403 | 0.954 | 0.535
Auto w/ iNat + WebCC† | 0.902 | 0.918 | 0.514 | 0.453 | 0.961 | 0.839
Auto w/ Expert GSV† | 0.929 | 0.929 | 0.740 | 0.686 | 0.967 | 0.854
Auto w/ Combined† | 0.931 | 0.921 | 0.662 | 0.744 | 0.971 | 0.907
Expert GSV‡ | 0.670 | 0.706 | 0.267 | 0.233 | 0.815 | 0.537
Table 5: Performance on remote sensing-based crop type classification for the
four major crop types in Thailand. *Baseline refers to classifying all samples
as the most common class, rice. †Ground references are automatically labeled
using CNNs trained on the specified datasets ($n=81,000$). ‡Ground references
contain only expert-labeled GSV images ($n=984$). Error is obtained over 10
runs with 10 different bootstrapped training sets.

#### Street-view crop type classification

Our experiments on the 5 datasets showed that weak supervision with online
images can successfully classify crop type in street-view images (Table 4).
WebCC ($n=2,754$) on its own achieved 82.9% overall accuracy, similar to
iNaturalist ($n=8,638$) at 83.2%. For comparison, a baseline that classifies
all samples as the most common class (rice) would achieve 70.3% accuracy.
Cassava, despite possessing visual features that distinguish it from the other
crops, had the lowest F1 score under both WebCC and iNaturalist (9.8% and
11.0%, respectively), followed by maize (40.3% and 41.3%) and sugarcane (53.6%
and 54.2%). The low performance is likely due to the low fraction of street-
view-like field images for iNaturalist (5% for cassava; Table A.2) and label
noise for WebCC, as Google Images has fewer images of cassava so it returns
images of related keywords (e.g., cassava root). When merged together,
iNaturalist + WebCC ($n=11,392$) overall accuracy improved to 93.8%. F1 scores
for cassava, maize, and sugarcane also substantially improved. This
complementarity could be due to iNaturalist having higher label accuracy but
more out-of-domain images (e.g. closeups), while WebCC has lower label
accuracy but considerably more street-view-like images.

The CNN trained on the expert-labeled dataset ($n=947$) achieved 93.9%
accuracy, surpassed only by Combined ($n=12,339)$, made up of all 3, which
achieved 95.9%. The comparable performance of iNaturalist + WebCC versus these
labelled datasets shows the weak supervision approach can help minimize the
need for expert labelling.

#### Role of sliding window classification and probability threshold

The MHP approach to remove low-confidence classifications (Section 2.3) proved
to increase the accuracy and F1 score across all models for thresholds set to
0.7 or above. The top performing model trained on the Combined dataset
achieved 81.6% accuracy when directly classifying the whole image, 90.6% with
sliding windows but no MHP threshold, and 95.9% with a 0.9 MHP threshold.
Improvements were similar for the other training sets.

We note that a higher MHP threshold $\tau$ leads to some images being dropped
from the inference set, because sometimes no sliding windows have softmax
probabilities exceeding $\tau$. Despite the slight loss of data, we observed
that F1 score increased monotonically with $\tau$ all the way up to
$\tau=0.95$ (Figure A.2).

###  3.3 Country-wide crop type map

#### Automated vs. expert-labeled GSV ground references

We found that training a random forest on large CNN-generated GSV ground
references ($n=81,000$) resulted in crop type maps that were significantly
more accurate than training on small expert-labeled ground references
($n=984$), despite the CNN-based labels being imperfect (Table 5). Training on
expert-labeled ground references achieved 67.0% overall accuracy, which is
actually lower than the 69.5% accuracy of a baseline model that classifies
everything as the most common class (rice). In comparison, the random forest
trained on 81,000 GSV samples whose crop types were predicted by the street-
view CNN trained on expert labels achieved an accuracy of 92.9%. This suggests
that the sample volume generated by leveraging GSV outweighs the noise that
predicted crop type labels can carry from misclassification.

#### Weakly supervised vs. expert-labeled automated GSV ground references

Among the 5 automatically-generated GSV ground reference datasets, the
datasets generated by CNNs trained on expert labels proved more accurate than
those generated by CNNs trained on weak labels from the web. Consistent with
the street-view classification results in Table 4, the lowest-performing
datasets were those created by the WebCC CNN and the iNaturalist CNN; the
remote sensing-based crop type classifier trained using their predictions as
labels achieved accuracies of 79.4% and 78.3%, respectively. The crop type
classifer trained on iNaturalist + WebCC CNN predictions performed better at
90.2% accuracy, but suffered from low cassava and maize F1 scores (51.4% and
45.3%). The two highest-performing crop type classifiers were trained on
outputs of the expert-labeled CNN and Combined CNN, achieving accuracies of
92.9% and 93.1% and higher F1 scores across all four crop types.

For the individual crop types, rice—the most abundant crop—was classified most
accurately, with F1 scores over 96% for the top 3 models. Sugarcane was also
classified accurately, with F1 scores ranging from 84% to 91% for the top 3
models. Maize and cassava were consistently the most difficult crop types to
classify across all models.

#### Sensitivity to training set size

Lastly, we investigated the relationship between the number of ground
references and the remote sensing-based crop type classification performance.
We found that, even at 81,000 ground references, performance had not yet
saturated as a function of training set size; the overall F1 score continued
to increase linearly as more street-view points were added (Figure 6).

Figure 6: Crop type map F1 score vs. ground reference dataset size. Computed
for the Auto w/ iNat + WebCC ground reference set from Table 5.

##  4 Discussion

We show that deep learning and street-view images can be combined to generate
thousands of geolocated crop type ground references at scale in smallholder
regions. These ground references, despite containing some label noise, can
then be used to create high-accuracy crop type maps in countries where no such
maps currently exist. In Thailand, 81,000 automated ground references led to a
more accurate crop type map than 1000 expert-labeled ground references, and
even at 81,000 performance had not yet saturated as a function of sample size.

To minimize the need for experts to manually label crop types in street-view
images, we explored using images from the web to weakly supervise a CNN to
classify crop type. We found that images from Google Images and iNaturalist,
although high in noise and often off-domain, can successfully supervise the
classification of street-view images. Furthermore, combining images from
different online sources improved performance. When creating the country-wide
crop type map in Thailand, we did observe that weakly supervised CNNs led to
lower accuracy crop type maps than CNNs trained on expert labels, suggesting
that the best solution may be to combine noisy images from the web with a
small number of expert-labeled street-view images.

One limitation of our work is that, to train a classifier to remove street-
view images containing small bushes, we manually labeled a set of field/not
field images. However, we point out that, unlike crop type labeling, field/not
field labeling does not require domain expertise. Another limitation is the
uncertain update frequency of street-view services and their continued uneven
distribution around the globe. Since more countries have street-view imagery
than crop type maps, we believe street-view datasets still have significant
value to add to global crop type mapping.

We release all datasets used for training and testing each model, including
over 81,000 crop type ground reference points, and the code to run this
pipeline in other regions. We also release the crop type map created for
Thailand’s 2022 wet growing season, which currently does not exist.

## References

  * Azzari and Lobell (2017) Azzari, G.; and Lobell, D. 2017.  Landsat-based classification in the cloud: An opportunity for a paradigm shift in land cover monitoring.  _Remote Sensing of Environment_ , 202: 64–74.
  * Becker-Reshef et al. (2023) Becker-Reshef, I.; Barker, B.; Whitcraft, A.; Oliva, P.; Mobley, K.; Justice, C.; and Sahajpal, R. 2023.  Crop Type Maps for Operational Global Agricultural Monitoring.  _Scientific Data_ , 10(1): 172.
  * Biljecki and Ito (2021) Biljecki, F.; and Ito, K. 2021.  Street view imagery in urban analytics and GIS: A review.  _Landscape and Urban Planning_ , 215: 104217.
  * Breiman (2001) Breiman, L. 2001.  Random forests.  _Machine learning_ , 45: 5–32.
  * Cantelaube and Carles (2014) Cantelaube, P.; and Carles, M. 2014.  Le registre parcellaire graphique: des données géographiques pour décrire la couverture du sol agricole.  _Le Cahier des Techniques de l’INRA_ , 58–64.
  * Carlini and Wagner (2017) Carlini, N.; and Wagner, D. 2017.  Towards evaluating the robustness of neural networks.  In _2017 ieee symposium on security and privacy (sp)_ , 39–57. Ieee.
  * Economic et al. (2014) Economic, F.; et al. 2014.  The State of Food and Agriculture (SOFA) 2014.
  * FAO (2021) FAO. 2021.  FAOSTAT [electronic resource].  https://www.fao.org/faostat.  Accessed: 2023-05-05.
  * Fisette et al. (2013) Fisette, T.; Rollin, P.; Aly, Z.; Campbell, L.; Daneshfar, B.; Filyer, P.; Smith, A.; Davidson, A.; Shang, J.; and Jarvis, I. 2013.  AAFC annual crop inventory.  In _2013 Second International Conference on Agro-Geoinformatics (Agro-Geoinformatics)_ , 270–274. IEEE.
  * Ghazaryan et al. (2018) Ghazaryan, G.; Dubovyk, O.; Löw, F.; Lavreniuk, M.; Kolotii, A.; Schellberg, J.; and Kussul, N. 2018.  A rule-based approach for crop identification using multi-temporal and multi-sensor phenological metrics.  _European Journal of Remote Sensing_ , 51(1): 511–524.
  * Gislason, Benediktsson, and Sveinsson (2006) Gislason, P. O.; Benediktsson, J. A.; and Sveinsson, J. R. 2006.  Random forests for land cover classification.  _Pattern recognition letters_ , 27(4): 294–300.
  * Haasnoot et al. (2013) Haasnoot, M.; Kwakkel, J. H.; Walker, W. E.; and Ter Maat, J. 2013.  Dynamic adaptive policy pathways: A method for crafting robust decisions for a deeply uncertain world.  _Global environmental change_ , 23(2): 485–498.
  * Han et al. (2012) Han, W.; Yang, Z.; Di, L.; and Mueller, R. 2012.  CropScape: A Web service based application for exploring and disseminating US conterminous geospatial cropland data products for decision support.  _Computers and Electronics in Agriculture_ , 84: 111–123.
  * Huerta-Ramos and Luštrik (2021) Huerta-Ramos, G.; and Luštrik, R. 2021.  Inat_Images: v1.1.
  * Inglada et al. (2015) Inglada, J.; Arias, M.; Tardy, B.; Hagolle, O.; Valero, S.; Morin, D.; Dedieu, G.; Sepulcre, G.; Bontemps, S.; Defourny, P.; et al. 2015.  Assessment of an operational system for crop type map production using high temporal and spatial resolution satellite optical imagery.  _Remote Sensing_ , 7(9): 12356–12379.
  * Jakubauskas et al. (2001) Jakubauskas, M. E.; Legates, D. R.; Kastens, J. H.; et al. 2001.  Harmonic analysis of time-series AVHRR NDVI data.  _Photogrammetric engineering and remote sensing_ , 67(4): 461–470.
  * Jin et al. (2019) Jin, Z.; Azzari, G.; You, C.; Di Tommaso, S.; Aston, S.; Burke, M.; and Lobell, D. B. 2019.  Smallholder maize area and yield mapping at national scales with Google Earth Engine.  _Remote Sensing of Environment_ , 228: 115–128.
  * Krause et al. (2016) Krause, J.; Sapp, B.; Howard, A.; Zhou, H.; Toshev, A.; Duerig, T.; Philbin, J.; and Fei-Fei, L. 2016.  The unreasonable effectiveness of noisy data for fine-grained recognition.  In _Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part III 14_ , 301–320. Springer.
  * Lee et al. (2022) Lee, J. Y.; Wang, S.; Figueroa, A. J.; Strey, R.; Lobell, D. B.; Naylor, R. L.; and Gorelick, S. M. 2022.  Mapping Sugarcane in Central India with Smartphone Crowdsourcing.  _Remote Sensing_ , 14(3): 703.
  * Lobell and Gourdji (2012) Lobell, D. B.; and Gourdji, S. M. 2012.  The influence of climate change on global crop productivity.  _Plant physiology_ , 160(4): 1686–1697.
  * Lowder, Skoet, and Raney (2016) Lowder, S. K.; Skoet, J.; and Raney, T. 2016.  The number, size, and distribution of farms, smallholder farms, and family farms worldwide.  _World development_ , 87: 16–29.
  * Ok, Akar, and Gungor (2012) Ok, A. O.; Akar, O.; and Gungor, O. 2012.  Evaluation of random forest method for agricultural crop classification.  _European Journal of Remote Sensing_ , 45(1): 421–432.
  * OpenStreetMap contributors (2017) OpenStreetMap contributors. 2017.  Planet dump retrieved from https://planet.osm.org .  https://www.openstreetmap.org.
  * Paliyam et al. (2021) Paliyam, M.; Nakalembe, C.; Liu, K.; Nyiawung, R.; and Kerner, H. 2021.  Street2sat: A machine learning pipeline for generating ground-truth geo-referenced labeled datasets from street-level images.  In _ICML 2021 Workshop on Tackling Climate Change with Machine Learning_.
  * Rezaei et al. (2021) Rezaei, E. E.; Ghazaryan, G.; Moradi, R.; Dubovyk, O.; and Siebert, S. 2021.  Crop harvested area, not yield, drives variability in crop production in Iran.  _Environmental Research Letters_ , 16(6): 064058.
  * Shumway, Stoffer, and Stoffer (2000) Shumway, R. H.; Stoffer, D. S.; and Stoffer, D. S. 2000.  _Time series analysis and its applications_ , volume 3.  Springer.
  * Szegedy et al. (2013) Szegedy, C.; Zaremba, W.; Sutskever, I.; Bruna, J.; Erhan, D.; Goodfellow, I.; and Fergus, R. 2013.  Intriguing properties of neural networks.  _arXiv preprint arXiv:1312.6199_.
  * Tseng et al. (2021) Tseng, G.; Zvonkov, I.; Nakalembe, C. L.; and Kerner, H. 2021.  Cropharvest: A global dataset for crop-type classification.  In _Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2)_.
  * Van Tricht et al. (2023) Van Tricht, K.; Degerickx, J.; Gilliams, S.; Zanaga, D.; Savinaud, M.; Battude, M.; Buguet de Chargère, R.; Dubreule, G.; Grosu, A.; Brombacher, J.; Pelgrum, H.; Lesiv, M.; Bayas, J. C. L.; Karanam, S.; Fritz, S.; Becker-Reshef, I.; Franch, B.; Bononad, B. M.; Cintas, J.; Boogaard, H.; Pratihast, A. K.; Kucera, L.; and Szantoi, Z. 2023.  ESA WorldCereal 10 m 2021 v100.
  * Wang, Azzari, and Lobell (2019) Wang, S.; Azzari, G.; and Lobell, D. B. 2019.  Crop type mapping without field-level labels: Random forest transfer and unsupervised clustering techniques.  _Remote sensing of environment_ , 222: 303–317.
  * Wang et al. (2020) Wang, S.; Di Tommaso, S.; Faulkner, J.; Friedel, T.; Kennepohl, A.; Strey, R.; and Lobell, D. B. 2020.  Mapping crop types in southeast India with smartphone crowdsourcing and deep learning.  _Remote Sensing_ , 12(18): 2957.
  * Wei et al. (2021) Wei, H.; Tao, L.; XIE, R.; and An, B. 2021.  Open-set Label Noise Can Improve Robustness Against Inherent Label Noise.  In Ranzato, M.; Beygelzimer, A.; Dauphin, Y.; Liang, P.; and Vaughan, J. W., eds., _Advances in Neural Information Processing Systems_ , volume 34, 7978–7992. Curran Associates, Inc.
  * Wu et al. (2021) Wu, F.; Wu, B.; Zhang, M.; Zeng, H.; and Tian, F. 2021.  Identification of crop type in crowdsourced road view photos with deep convolutional neural network.  _Sensors_ , 21(4): 1165.
  * Yan and Ryu (2021) Yan, Y.; and Ryu, Y. 2021.  Exploring Google Street View with deep learning for crop type mapping.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 171: 278–296.

##  Appendix A Appendix

###  A.1 Generating equidistant points from OSM

Line segments are created for each OSM way, by creating an edge between OSM
node pairs using linear interpolation. Then starting with the point at the
base of the line, latitude and longitude points are generated at 10 meter
steps recursively over the distance between the node pairs as shown below:

$\text{lat}_{\text{new}}=\text{lat}_{\text{prev}}+\frac{\text{distance}}{10\cdot(\text{lat}_{\text{new}}-\text{lat}_{\text{prev}})}$
(1)

$\text{lon}_{\text{new}}=\text{lon}_{\text{prev}}+\frac{\text{distance}}{10\cdot(\text{lon}_{\text{new}}-\text{lon}_{\text{prev}})}$
(2)

###  A.2 Pulling GSV street-level imagery

Google Street View (GSV) is used as the dataset to gather images at the
desired geocoordinates in fields. To query the GSV API, the arguments (lat,
lon) and camera heading ($\theta$) are used to download images at the desired
location with the camera facing towards the field. GSV charges $7 for 1000
images, which is why a significant filtering on the set of points is performed
beforehand using land cover maps. The street bearing is necessary for two
purposes (1) to pull GSV images with the camera facing the direction of the
field, that is, +90 degrees and -90 degrees from the street bearing and (2) to
calculate the point coordinate of labelled fields that will be fed to the
remote sensing model.

To achieve this the road bearing for each filtered street point is
interpolated using Harversine formulas - which account for earth’s ellipsoidal
shape. The formulas calculate the direction $\theta$ from a start point
$(\text{lat}_{1},\text{lon}_{1})$ to its adjacent point
$(\text{lat}_{2},\text{lon}_{2})$ 10m away. Latitudes and longitudes are in
radians.

$\Delta lat=lat_{2}-lat_{1}$ (3)
$y=\operatorname{sin}(lat_{2}-lat_{1})\cdot\operatorname{cos}(lon_{2})$ (4)
$x=\operatorname{cos}(lon_{1})\cdot\operatorname{sin}(lon_{2})-\operatorname{sin}(lon_{1})\cdot\operatorname{cos}(lon_{2})\cdot\operatorname{cos}(\Delta
lat)$ (5) $\theta=\operatorname{atan}\tfrac{y}{x}$ (6)

Once the street bearing is calculated and converted back to degrees, the
bearing for the camera to face the two fields on either side of the street are
$\theta$ \+ 90 and $\theta$ \- 90.

Next we calculate the field point D meters away, perpendicular to the street
bearing, where a GSV image was taken. This field point will be labelled in
section 4.2 and used as ground reference in training data for remote sensing
in section 4.3. To achieve this the direct haversine formula is used.
$\textit{lat}_{s}$ and $\textit{lon}_{s}$ are street point coordinates to
calculate, $\textit{lat}_{f}$ and $\textit{lon}_{f}$ are field point
coordinates, D is the distance away from the street, $\theta$ is the
direction, R is the mean radius of the earth, and $\textit{A}_{D}$ =
$\frac{D}{R}$.

$lon_{f}=\operatorname{asin}(\sin(lon_{s})\cdot\cos(A_{D})+\cos(lon_{s})\cdot\sin(A_{d})\cdot\cos(\theta))$
(7)

$lat_{f}=lat_{s}+\operatorname{atan}(\tfrac{\sin(\theta)\cdot\operatorname{sin}(A_{d})\cdot\operatorname{cos}(lon_{s})}{\operatorname{cos}(A_{D})-\operatorname{sin(lon_{s})}\cdot\operatorname{sin}(lon_{f})})$
(8)

The $(lon_{f},lat_{f})$ field geocoordinates are assigned to their respective
image to be labelled, to serve together as a ground reference point to train
the remote sensing models.

Distance (m) | Furthest FB% | Closest FB%
---|---|---
<10 | 0 | 8
10 | 0 | 83
20 | 0 | 9
30 | 2 | 0
40 | 21 | 0
50 | 41 | 0
>50 | 36 | 0
Table A.1: Showing the distance of furthest field boundary from to the street, and the distance of the closest boundary from the street. The data was collected by inspecting random samples across Thailand for 200 fields. FB stands for field boundary. The experiment showed placing the ground reference 30 m perpendicular from where the GSV image was taken is reasonable with only 2% of outliers.  Figure A.1: OSM Points in black and equidistant generated points in red. OSM points gathered from OSM ways are random and sparse. The red points are generated at 10m steps in order to query GSV and ensure higher recall of field points and street-view images.  Figure A.2: Solid lines show the impact of increasing the Mode of High Probs (MHP) threshold on the test set F1 score for the different models. Circles represent the sampled points. Dotted lines show the impact of increasing the MHP threshold on the percent of GSV field images being labelled that are dropped from the ground referenced set because they contain no softmax probabilities over the threshold.  Figure A.3: Example of the raw GCVI time series for the four major crops in Thailand (cassava, maize, sugarcane and rice). Dataset | Cassava # (SV%) | Maize # (SV%) | Rice # (SV%) | Sugarcane # (SV%) | Total # (SV%)
---|---|---|---|---|---
iNaturalist | 3662 (5%) | 1698 (19%) | 1396 (34%) | 1882 (20%) | 8646 (16%)
WebCC | 584 (77%) | 832 (65%) | 659 (48%) | 679 (66%) | 2754 (71%)
Table A.2: Distribution of scraped online datasets by crop type. SV% stands
for percentage of street view like images in the online datasets. CC stands
for creative commons.


#  Low-rank quantics tensor train representations of Feynman diagrams for
multiorbital electron-phonon models

Hirone Ishida  Department of Physics, Saitama University, Saitama 338-8570,
Japan  Natsuki Okada  Department of Physics, Saitama University, Saitama
338-8570, Japan  Shintaro Hoshino  Department of Physics, Saitama University,
Saitama 338-8570, Japan  Hiroshi Shinaoka  Department of Physics, Saitama
University, Saitama 338-8570, Japan

###### Abstract

Feynman diagrams are an essential tool for simulating strongly correlated
electron systems. However, stochastic quantum Monte Carlo (QMC) sampling
suffers from the sign problem, e.g., when solving a multiorbital quantum
impurity model. Recently, two approaches have been proposed for efficient
numerical treatment of Feynman diagrams: Tensor Cross Interpolation (TCI) for
replacing the stochastic sampling and the Quantics Tensor Train (QTT)
representation for compressing space-time dependence. Combining these
approaches, we find low-rank structures in weak-coupling Feynman diagrams for
a multiorbital electron-phonon model and demonstrate their efficient numerical
integrations with exponential resolution in time and exponential convergence
of error with respect to computational cost.

In diverse fields ranging from condensed matter physics to high-energy
physics, calculations based on Feynman diagrams in quantum field theory have
become indispensable [1]. Numerical perturbative calculations often involve
numerical integrations and summations over virtual times and internal degrees
of freedom of these diagrams. However, when dealing with multiorbital models
or higher-order diagrams, straightforward numerical integration faces the
notorious “curse of dimensionality,” leading to exorbitant computational
times. Consequently, methods to circumvent this curse in quantum field theory
calculations are in high demand.

Various approaches have been explored to address this challenge. For instance,
quantum Monte Carlo (QMC) methods [2, 3, 4, 5], grounded in importance
sampling in multidimensional spaces, are widely employed. Yet, in complex
systems, these methods suffer from the negative sign problem, causing
computational times to explode exponentially [6, 7, 8, 9, 10]. This
occasionally hinders accessing parameter regimes with physically interesting
phenomena.

A recent breakthrough has been the introduction of Tensor Cross Interpolation
(TCI), which employs tensor networks to learn a function and perform its
integration in multidimensional spaces [11]. The TCI approach is based on low-
rank structures hidden in a function. This method has demonstrated its prowess
by calculating physical quantities with unprecedented precision in various
models, outperforming the Monte Carlo methods in some cases, e.g., in solving
single-orbital quantum impurity models [11, 12]. They observed that the
integration error in TCI converges more rapidly than in QMC sampling,
following a rate of $O(1/n_{\mathrm{s}}^{a})$, with $a>1/2$, where
$n_{\mathrm{s}}$ denotes the number of integrand evaluations. Yet, determining
the factors influencing this convergence rate remains unresolved. Furthermore,
extending the TCI method to multiorbital systems and electron-phonon coupling
remains to be developed, which is essential for its practical application in
materials calculations.

Figure 1:  (a) Second-order skeleton diagrams for the electron self-energy
$\Sigma^{(2)}$ and phonon self-energy $\Pi^{(2)}$, as well as tensor train
(TT) representations of their integrands. (b) Bond dimensions $\chi_{\ell}$
for TCI tolerance $\epsilon=10^{-3}$ for typical parameters set above/below
the superconducting temperature $T_{\mathrm{c}}\simeq 0.2$ and with/without an
external field, whose strength is denoted by $\alpha$. We employed a virtual
$2^{{\mbox{$\mathcal{R}$}}}\times 2^{{\mbox{$\mathcal{R}$}}}\times
2^{{\mbox{$\mathcal{R}$}}}$ grid for $\tau$, $\tau^{\prime}$ and
$\tau^{\prime\prime}$ (${{\mbox{$\mathcal{R}$}}}=20$,
$2^{{\mbox{$\mathcal{R}$}}}\simeq 10^{6}$). (c) Interpolation error estimate
as a function of bond dimension $\chi$. The error is normalized by an estimate
of the absolute maximum of the function (integrand).

In a complementary vein, a novel discretization technique for space-time based
on the quantics tensor train (QTT) has been proposed for solving Navier-Stokes
equations for turbulent flows [13, 14], the Vlasov-Poisson equations for
collisionless plasmas [15], and solving diagrammatic equations in quantum
field theories [16]. The QTT method leverages the low-entanglement structures
between exponentially different length scales, allowing for an exponential
enhancement in spatiotemporal resolution relative to data volume and
computational effort. In Ref. [16], the separation in length scales has been
showcased in various multidimensional spacetime datasets of correlation
functions in quantum field theories, including momentum dependence of one-
particle Green’s function, three-frequency dependence of a vertex function at
the two-particle level and so on.

Furthermore, the amalgamation of TCI and QTT has given birth to the quantics
tensor cross interpolation (QTCI) technologies [17]. They illustrated its
potential with an application from condensed matter physics: the computation
of Brillouin zone integrals [17]. Applications of QTCI to more complex cases,
e.g., integration over both continuous and discrete variables, are promising
and remain to be explored.

The primary objective of this study is to address the following questions: (i)
Do low-rank tensor train (TT) structures exist within self-energy Feynman
diagrams for multiorbital systems? (ii) Is the TCI robust for interpolating
complex integrands with continuous and discrete degrees of freedom? (iii) Can
TCI accommodate phonon contributions, which, despite their physical relevance,
occasionally worsens the sign problem in QMC simulations? (iv) Does employing
quantics facilitate achieving an enhanced convergence rate, specifically
exponential convergence? To investigate these questions, as a first step, we
analyze weak-coupling self-energy Feynman diagrams in the imaginary-time
formalism for a prototype multiorbital electron-phonon model, originally
proposed for exploring superconducting states in fullerides [18, 19, 20, 21,
22].

Figure 1 summarizes our main results obtained with and without an external
field. Figure 1(a) illustrates the TT structures proposed in this study for
the electron and phonon self-energies, where imaginary-time variables are
discretized by the quantics representation. Discrete internal and external
variables, such as spin-orbital indices $\gamma$ and vibrational modes of
phonon $\eta$ are encoded in the same TT as well. As shown in Fig. 1(b), the
bond dimensions of the TT representations obtained by TCI are much smaller
than the worst case of incompressible functions, establishing the existence of
low-rank structures. As demonstrated in Fig. 1(c), the interpolation error
estimate in the TCI vanishes exponentially regardless of the superconducting
transition and the presence of the external field, indicating the advantage of
the QTT representation.

The Hamiltonian of the model reads

$\displaystyle\mathscr{H}$
$\displaystyle=\sum_{ij}\sum_{\gamma\gamma^{\prime}\sigma}\left(t^{\gamma\gamma^{\prime}}_{ij}-\mu\delta_{ij}\delta_{\gamma\gamma^{\prime}}\right)c^{\dagger}_{i\gamma\sigma}c_{j\gamma^{\prime}\sigma}+\sum_{i\eta}\omega_{\eta}a^{\dagger}_{i\eta}a_{i\eta}$
$\displaystyle\quad+\sum_{i\eta}I_{\eta}:T_{i\eta}T_{i\eta}:+\sum_{i\eta}g_{\eta}\phi_{i\eta}T_{i\eta},$
(1)

where $\gamma=x,y,z$ represent the $t_{1u}$ electronic orbitals,
$\eta=0,1,3,4,6,8$ correspond to the six vibrational modes of a fullerene
molecule with energies $\omega_{\eta}$, and $c^{\dagger}$ and $c$ are the
creation and annihilation operators for electrons, respectively, while
$a^{\dagger}$ and $a$ are the creation and annihilation operators for phonons.
The $\lambda^{\eta}$ are Gell-Mann matrices,
$T_{i\eta}=\sum_{\gamma\gamma^{\prime}\sigma}c^{\dagger}_{i\gamma\sigma}\lambda^{\eta}_{\gamma\gamma^{\prime}}c_{i\gamma^{\prime}\sigma}$,
$\phi_{i\eta}=a_{i\eta}+a^{\dagger}_{i\eta}$, $I_{\eta}=(3/4)U-J$ ($\eta=0$),
$=J/2$ (otherwise) . The first term represents inter-site hopping $t_{ij}$ and
the chemical potential $\mu$. The non-interacting part of the Hamiltonian is
given by the semicircular non-interacting density of states of width 1. The
second term gives the single-particle phonon energies. The third term
represents instantaneous local Columb interactions, where the ($:$) symbol
represents the normal ordering. The fourth term denotes the electron-phonon
couplings. For more details on the model, we refer the reader to Ref. 22.

We take $\omega_{\eta}=\omega_{0}=0.15$ and
$g_{\eta}=g_{0}=\sqrt{3\lambda_{0}\omega_{0}/4}$ with $\lambda_{0}=0.15$. We
set $U=J=\mu=0$ and consider only the electron-phonon interaction to
investigate low-rank structures hidden in Feynmann diagrams functions.
Numerical treatment of the instantaneous part will be more straightforward.

Figure 2:  (a) $T$ dependence of the absolute maximum of the gap function with
no external field ($\alpha=0$). (b) $\alpha$ dependence of the absolute
maximum of the gap function at $T=0.016$. (c) $\Sigma^{(1)}$ and $\Pi^{(1)}$
computed at $T=0.05,~{}0.016$ with $\alpha=0$. (d) $\Sigma^{(2)}$ and
$\Pi^{(2)}$ computed at $T=0.05,~{}0.016$ with $\alpha=0$.

We solve the model within the dynamical-mean field theory by assuming local
self-energies. The local Green’s function in the Nambu basis reads

$\displaystyle\check{G}_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}(\tau)$
$\displaystyle=-\expectationvalue{\mathcal{T}\Psi(\tau)\Psi^{\dagger}(0)}=\left(\begin{smallmatrix}G_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}(\tau)&F_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}(\tau)\\\
\bar{F}_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}(\tau)&\bar{G}_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}(\tau)\end{smallmatrix}\right),$
$\displaystyle D_{\eta\eta^{\prime}}(\tau)$
$\displaystyle=-\expectationvalue{\mathcal{T}\phi_{\eta}(\tau)\phi_{\eta^{\prime}}(0)},$

where
$\Psi(\tau)=(c_{x\uparrow}(\tau),c_{x\downarrow}(\tau),\cdots,c_{z\downarrow}(\tau),c^{\dagger}_{x\uparrow}(\tau),\cdots,c^{\dagger}_{z\downarrow}(\tau))^{\\!\top}$,
$c_{\gamma\sigma}(\tau)$/$c_{\gamma\sigma}^{\dagger}(\tau)$ represents an
annihilation/creation operator for orbital $\gamma$ and spin $\sigma$ in the
Heisenberg picture, $\mathcal{T}$ indicates time odering.

Note that matrices with a checkmark $\check{}$, such as $\check{A}$, are
defined in the Nambu space and have a shape of $12\times 12$.

A weak-coupling expansion yields the following connected skeleton diagrams for
the electron and phonon self-energies [Fig. 1(a)]:

$\displaystyle\check{\Sigma}^{(1)}(\tau)=-\sum_{\eta_{1}\eta_{2}}g_{0}^{2}D_{\eta_{1}\eta_{2}}(\tau)\check{\lambda}^{\eta_{1}}\check{G}(\tau)\check{\lambda}^{\eta_{2}},$
(2)
$\displaystyle\check{\Sigma}^{(2)}(\tau)=\int_{[0,\beta]^{2}}d\tau^{\prime}d\tau^{\prime\prime}\sum_{\eta_{1}\eta_{2}\eta_{3}\eta_{4}}g_{0}^{4}D_{\eta_{1}\eta_{3}}(-\tau^{\prime\prime})D_{\eta_{2}\eta_{4}}(\tau^{\prime}-\tau)$
$\displaystyle~{}~{}~{}\times\check{\lambda}^{\eta_{4}}\check{G}(\tau-\tau^{\prime\prime})\check{\lambda}^{\eta_{3}}\check{G}(\tau^{\prime\prime}-\tau^{\prime})\check{\lambda}^{\eta_{2}}\check{G}(\tau^{\prime})\check{\lambda}^{\eta_{1}},$
(3)
$\displaystyle\Pi^{(1)}_{\eta\eta^{\prime}}(\tau)=\dfrac{1}{2}g_{0}^{2}\Tr\left[\check{G}(-\tau)\check{\lambda}^{\eta^{\prime}}\check{G}(\tau)\check{\lambda}^{\eta}\right],$
(4)
$\displaystyle\Pi^{(2)}_{\eta\eta^{\prime}}(\tau)=-\dfrac{1}{2}\int_{[0,\beta]^{2}}d\tau^{\prime}d\tau^{\prime\prime}\sum_{\eta_{2}\eta_{4}}g_{0}^{4}D_{\eta_{2}\eta_{4}}(\tau^{\prime}-\tau^{\prime\prime})\times$
$\displaystyle~{}{\rm
Tr}\left[\check{G}(-\tau^{\prime\prime})\check{\lambda}^{\eta_{4}}\check{G}(\tau^{\prime\prime}-\tau)\check{\lambda}^{\eta}\check{G}(\tau-\tau^{\prime})\check{\lambda}^{\eta_{2}}\check{G}(\tau^{\prime})\check{\lambda}^{\eta^{\prime}}\right].,$
(5)

where $\beta=1/T$ denotes inverse temperature. The Gell-Mann matrix in the
Nambu space reads
$\check{\lambda}^{\eta}=\mathcal{N}(\lambda^{\eta}\otimes\mathbf{1}_{2\times
2})$, where $\mathbf{1}_{2\times 2}$ is the unit matrix in spin space
($2\times 2$) and
$\check{M}=\mathcal{N}(M)\equiv\left(\begin{smallmatrix}M&0\\\
0&-M^{\\!\top}\end{smallmatrix}\right)$. The electron’s self-energy, denoted
by $\check{\Sigma}$, is defined as a $12\times 12$ matrix in the form of

$\displaystyle\check{\Sigma}(\tau)=\begin{pmatrix}\Sigma(\tau)&\Delta(\tau)\\\
\Delta^{\dagger}(\tau)&-\Sigma^{\\!\top}(-\tau)\end{pmatrix}\equiv\begin{pmatrix}\left(\check{\Sigma}\right)_{11}&\left(\check{\Sigma}\right)_{12}\\\
\left(\check{\Sigma}\right)_{21}&\left(\check{\Sigma}\right)_{22}\\\
\end{pmatrix},$ (6)

where $\Delta(\tau)$, referred to as the gap function, is a $6\times 6$
matrix. The phonon Green’s function $D$ and the self-energy $\Pi$ are $6\times
6$ matrix.

The self-consistent equations read

$\displaystyle\\!\\!\\!\check{G}({\rm i}\omega_{n})\\!$
$\displaystyle=\\!\\!\int\\!d\varepsilon\rho(\varepsilon)\\!\left[{\rm
i}\omega_{n}\mathbf{1}_{12\times
12}\\!+\\!(\mu-\varepsilon)\check{\mathbf{1}}\\!-\\!\check{\Sigma}({\rm
i}\omega_{n})\\!-\\!\alpha\mathscr{H}_{\mathrm{ex}}\right]^{-1},$
$\displaystyle D({\rm i}\nu_{m})$ $\displaystyle=D^{(0)}({\rm
i}\nu_{m})[\mathbf{1}_{6\times 6}-D^{(0)}({\rm i}\nu_{m})\Pi({\rm
i}\nu_{m})]^{-1}$

where $\omega_{n}$ and $\nu_{m}$ are Matsubara frequencies for fermions and
bosons, $\mathbf{1}_{N\times N}$ is the $N\times N$ unit matrix,
$\check{\mathbf{1}}\equiv\mathcal{N}(\mathbf{1}_{6\times 6})$ and
$\rho(\varepsilon)=(8/\pi)\sqrt{(1/2)^{2}-\varepsilon^{2}}$ ,
$D_{\eta\eta^{\prime}}^{(0)}({\rm
i}\nu_{m})=2\delta_{\eta\eta^{\prime}}\omega_{0}/(({\rm
i}\nu_{m})^{2}-\omega_{0}^{2})$ is the non-interacting phonon Green’s
function.

We can apply an external field with $\alpha\neq 0$. We apply the external
field to lower the symmetry to examine cases where the symmetry is reduced.
Here, we chose
$\mathscr{H}_{\mathrm{ex}}=\mathcal{N}(l_{z}\otimes\mathbf{1}_{2\times 2})$ as
a test field. Here, $l_{z}=\left(\begin{smallmatrix}0&-{\rm i}&0\\\ {\rm
i}&0&0\\\ 0&0&0\end{smallmatrix}\right)$ represents the orbital magnetic
field. $\alpha$ denotes the magnitude of the external field.

We perform self-consistent calculations with $\Sigma^{(1)}$ and $\Pi^{(1)}$
using the sparse modeling methods based on the intermediate representation
(IR) [23, 24, 25] to describe a superconducting state. The resultant Green’s
functions are stored in IR and can be evaluated at any $\tau$. This first-
order solution is used for the evaluation of the more complex second-order
diagrams involving imaginary time integrals.

Quantics tensor train We start with introducing the quantics representation of
a univariate function, denoted by $f(x):x\in[0,1)\rightarrow\mathbb{C}$. The
function can be “tensorized” by binary notation [26, 27, 16]. In practice, we
construct an equidistant grid on the interval $[0,1)$:
$x=(0.\sigma_{1}\sigma_{2}\dots\sigma_{{{\mbox{$\mathcal{R}$}}}})_{2}=\sum_{r=1}^{{{\mbox{$\mathcal{R}$}}}}\sigma_{r}2^{-r}$,
where $\sigma_{r}\in\\{0,1\\}$ corresponds to the scale $2^{-r}$. We then
decompose the $\mathcal{R}$-way tensor,
$F_{\sigma_{1},\dots,\sigma_{{{\mbox{$\mathcal{R}$}}}}}\equiv
f((0.\sigma_{1}\cdots\sigma_{{\mbox{$\mathcal{R}$}}})_{2})$, as

$\displaystyle
F_{\sigma_{1},\sigma_{2},\cdots,\sigma_{{{{\mbox{$\mathcal{R}$}}}}}}\approx\sum_{\alpha_{1}}^{\chi_{1}}\cdots\sum_{\alpha_{{{{\mbox{$\mathcal{R}$}}}}-1}}^{\chi_{{{{\mbox{$\mathcal{R}$}}}}-1}}F^{(1)}_{\sigma_{1}\alpha_{1}\alpha_{2}}\cdots
F^{({{\mbox{$\mathcal{R}$}}})}_{\sigma_{{{{\mbox{$\mathcal{R}$}}}}}\alpha_{{{{\mbox{$\mathcal{R}$}}}}-1}},$
(7) $\displaystyle\hskip
20.00003pt\raisebox{-0.5pt}{\includegraphics[width=173.44534pt]{schematic_QTT.pdf}}~{},$
(8)

where $F_{\sigma_{1},\dots,\sigma_{{{\mbox{$\mathcal{R}$}}}}}$ has a total of
$2^{{\mbox{$\mathcal{R}$}}}$ elements, and $F^{(\ell)}$ is a three-way tensor
of size $2\times\chi_{\ell-1}\times\chi_{\ell}$. The sizes of the auxiliary
indices (virtual bonds) $\alpha_{\ell}$, $\chi_{\ell}$, are referred to as
bond dimensions. The bond dimension of the TT is defined as
$\chi\equiv\mathrm{max}_{\ell}~{}\chi_{\ell}$. The integral of $f(x)$ can be
evaluated accurately by Riemann summation as

$\displaystyle\raisebox{-0.5pt}{ \includegraphics[width=433.62pt]{TT_int.pdf}
}~{},$

where  denotes the vector $(1/2,1/2)$. In the above equation, the first
$\simeq$ denotes the approximation by the Riemann sum, where the
discretization error can be negligibly and exponentially small as
$O(2^{-2{{\mbox{$\mathcal{R}$}}}})$. The second $\simeq$ corresponds to the TT
approximation. The computational time of the contraction scales linearly with
$\mathcal{R}$ as $O(\chi^{2}{{\mbox{$\mathcal{R}$}}})$.

The exact representation of an incompressible function, for example, a random
function, requires $\chi_{\ell}$ to grow toward the center of the TT
exponentially, offering no memory advantage over
$F_{x_{1},\dots,x_{{{\mbox{$\mathcal{R}$}}}}}$. In the presence of a low-rank
structure in $F$, i.e., the separation between different length scales [16],
the bond dimension can be significantly reduced; as a result, the total number
of elements in the TT grows only linearly with $\mathcal{R}$ as
$\mathcal{O}(\chi^{2}{{\mbox{$\mathcal{R}$}}})$. The imaginary-time dependence
of propagators self-energies is expected to have such low-rank structures [16,
28].

TT representations for integrands The integrands of self-energies (3) and (5)
depend on multiple discrete and continuous variables. For example, the
integrand of the electron self-energy (3) depends on discrete variables
$\gamma,\gamma^{\prime},\eta_{1},\eta_{2},\eta_{3},\eta_{4}$, and continuous
variables $\tau,\tau^{\prime},\tau^{\prime\prime}$. We apply the binary coding
to the continuous variables as
$\tau/\beta=(0.\tau_{1}\cdots\tau_{{{\mbox{$\mathcal{R}$}}}})_{2}$,
$\tau^{\prime}/\beta=(0.\tau^{\prime}_{1}\cdots\tau^{\prime}_{{{\mbox{$\mathcal{R}$}}}})_{2}$,
$\tau^{\prime\prime}/\beta=(0.\tau^{\prime\prime}_{1}\cdots\tau^{\prime\prime}_{{{\mbox{$\mathcal{R}$}}}})_{2}$.
Figure 1(a) illustrates the TTs for the integrands. The lengths of the TTs are
${{\mbox{$\mathcal{L}$}}}=3{{\mbox{$\mathcal{R}$}}}+6$ and
${{\mbox{$\mathcal{L}$}}}=3{{\mbox{$\mathcal{R}$}}}+4$, respectively. We
assign $\tau_{r}$, $\tau^{\prime}_{r}$, $\tau^{\prime\prime}_{r}$, at the same
length scale to the same tensor due to the strong entanglement between these
degrees of freedom [16, 28]. The indices for the discrete variables are
grouped on the left side of the TTs.

Tensor cross interpolation (TCI) We construct a TT for the integrand by TCI,
which is an adaptive learning algorithm that constructs a low-rank TT from a
small subset of adaptively chosen elements of a tensor $T$ (interpolation
points). Specifically, it optimizes $\chi_{\ell}$ so that a global error
estimate satisfies
$\|T-\tilde{T}\|_{\mathrm{max}}/\|T\|_{\mathrm{max}}<\epsilon$, where
$\epsilon$ is a tolerance, $\tilde{T}$ is a TT and $\|\cdots\|_{\mathrm{max}}$
denotes the maximum norm. Generally, the number of interpolation points
required for constructing a TT with length $\mathcal{L}$ and bond dimension
$\chi$ scales linearly with $\mathcal{R}$ as $O(\chi^{2}\mathcal{L})$, leading
to an exponential speed-up over the conventional tensor decomposition based on
singular value decomposition. For details, please refer to Ref. [29].

For constructing a TT for the integrand, we evaluate the RHS of Eq. (3) and
Eq. (5) for indices chosen by TCI. This can be performed straightforwardly if
the $\tau$ dependence of the Green’s functions is expanded by IR. The function
evaluations consume the majority of the computational time of the TT
construction. We explored the configuration space by combining 2-site (local)
updates [29] and global moves [30]. We found that the use of the global moves
is necessary for a stable convergence of the TCI: The 2-site updates suffer
from an ergodicity problem [29].

Integration over internal variables Once the TT representation is constructed,
summation and integration over the internal variables can be performed
efficiently, as explained above. For instance, the TT representation of the
electron self-energy (3) can be constructed as 111 We must include the factor
$\beta^{2}$ due to changing variables: $x^{\prime}=\tau^{\prime}/\beta$ and
$x^{\prime\prime}=\tau^{\prime\prime}/\beta$.

$\displaystyle\raisebox{-0.5pt}{
\includegraphics[width=368.57964pt]{TT.pdf}},$

where the tensor with no leg (red) can be absorbed into either its left (for
$\gamma^{\prime}$) or right (for $\tau_{1}$) tensor. The resultant TT can be
readily evaluated at any point $(\gamma,\gamma^{\prime},\tau)$ with
exponentially high resolution for $\tau$.

We first discuss the $\alpha$-$T$ phase diagram. Figure 2(a) plots the $T$
dependence of the order parameter of the superconducting transition. The
system exhibits a second-order phase transition into a superconducting phase
around $T=0.02$, being consistent with the results in the previous study [22].
Figure 2(b) shows the $\alpha$ dependence of the order parameter at $T=0.016$.
The order parameter vanishes around $\alpha=0.08$, signaling a transition to
the normal phase.

Figure 3:  Integrands of (a) $\Sigma^{(2)}_{x\uparrow,x\uparrow}$, (b)
$\Delta_{x\uparrow,x\downarrow}^{(2)}$, and (c) $\Pi^{(2)}_{11}$ at
$\tau=\beta/2$ and $T=0.016$. Panels (d)–(f) show interpolation errors with
$\epsilon=10^{-3}$ for the integrands in (a)–(c), respectively.

Figure 2(c) shows typical self-energies $\Sigma^{(1)}$ and $\Pi^{(1)}$
computed at $T=0.05$ ($>T_{\mathrm{c}}$) and $T=0.016$ ($<T_{\mathrm{c}}$) for
$\alpha=0$. Let us discuss the structures of the self-energy matrices in the
flavor space. At $\alpha=0$ and $T>T_{\mathrm{c}}$, the normal part of the
electron self-energy $\Sigma^{(1)}(\tau)$ and the bosonic self-energy
$\Pi^{(1)}(\tau)$ are diagonal, while the anomalous part of the electron self-
energy $\Delta^{(1)}(\tau)=0$. At $\alpha=0$ and for $T<T_{\mathrm{c}}$, the
spontaneous symmetry breaking gives rise to nonzero superdiagonal and
subdiagonal elements in $\Delta^{(1)}(\tau)$. Introducing $\alpha>0$ further
reduces the symmetry. We observed that nonzero components appear in
$(\gamma\sigma,\gamma^{\prime}\sigma^{\prime})=(x\\!\\!\uparrow,y\\!\\!\uparrow),(x\\!\\!\downarrow,y\\!\\!\downarrow)$
for $\Sigma$,
$(\gamma\sigma,\gamma^{\prime}\sigma^{\prime})=(x\\!\\!\uparrow,y\\!\\!\downarrow),(x\\!\\!\downarrow,y\\!\\!\uparrow)$
for $\Delta$, and $(\eta,\eta^{\prime})=(1,3),(4,6),(0,8)$ for $\Pi$,
respectively. This refers to the upper triangular part, and nonzero components
also appear symmetrically in the lower triangular part. The appearance of
these nonzero components is consistent with what is expected from the symmetry
of the external field. The second-order self-energy matrices have precisely
the same structures of nonzero elements as those of the first-order self-
energies.

The left panels of Fig. 3 show typical components of the integrands for
$\Sigma^{(2)}$, $\Delta^{(2)}$, and $\Pi^{(2)}$ computed at $\alpha=0$ and
$T=0.016$. These integrands have several lines of discontinuities, originating
from the discontinuity of the one-particle propagators at $\tau=n\beta$
($n\in\mathbb{Z}$).

We now construct TCIs of these functions for typical parameter sets of $T$ and
$\alpha$. The tolerance is set to $\epsilon=10^{-3}$. The TCI construction
takes a few minutes with one Apple M2 CPU core for $\chi\simeq 100$.

The results are shown in Fig. 1(b). We numerically confirmed that the
constructed $\tau$-dependent $12\times 12$ matrix $\check{\Sigma}$ satisfies
the expected relation between the blocks:
$\left(\check{\Sigma}\right)_{11}(\tau)=-\left(\check{\Sigma}\right)_{22}^{\\!\top}(-\tau)$
and
$\left(\check{\Sigma}\right)_{12}(\tau)=\left(\check{\Sigma}\right)_{21}^{\dagger}(\tau)$
. This indicates the stability of the TCI.

Figure 1(b) shows the bond dimensions $\chi_{\ell}$ computed for typical
parameter sets. The bond dimensions are significantly smaller than the worst
(incompressible) case and strongly depend on $\ell$. For the electron self-
energy, $\chi_{\ell}$ takes large values at $\ell=4$ (between $\eta_{2}$ and
$\eta_{3}$) and at $\ell\simeq 8$. The bond $\ell=8$ connects the time scales
of $\tau/\beta=2^{-2}$ ($r=2$) and $2^{-3}$ ($r=3$). For shorter time scales,
$\chi_{\ell}$ monotonically decreases, allowing us to increase the time
resolution exponentially. We observed a similar trend for the phonon self-
energy: $\chi_{\ell}$ takes large values around $\ell=6$, corresponding to
between $r=2$ and $r=3$.

We move on to a more detailed discussion on the compactness of the quantics
representation. Figure 1(c) shows the $\chi$ dependence of an estimate of the
interpolation error. The interpolation error for all the parameter sets
eventually decreases exponentially with $\chi$. Let us first focus on the
cases with no external field ($\alpha=0$) in more depth. Comparing the results
for $T=0.02~{}(>T_{\mathrm{c}})$ and $T=0.016~{}(<T_{\mathrm{c}})$, one can
see that the symmetry breaking roughly doubles $\chi$ required to reach
$\epsilon=10^{-2}$ from $\simeq 35$ to $\simeq 55$. This can be attributed to
the emergence of the nonzero components in $\check{\Sigma}$.

At $T=0.016$, increasing $\alpha=0$ to $\alpha=0.05$ in the superconducting
phase does not change the behavior of the interpolation error above
$\epsilon=10^{-2}$ qualitatively but leads to the emergence of a more slowly
decaying tail below $\epsilon=10^{-2}$. This tail may originate from the
requirement of interpolating the small off-diagonal components in
$\check{\Sigma}$ induced by the external field. Further increasing $\alpha$
beyond the transition point, $\chi$ becomes substantially smaller with the
disappearance of the spontaneous symmetry breaking.

The right panels of Fig. 3 show the interpolation error in the reconstructed
self-energy for a typical parameter set, $T=0.016$ and $\alpha=0$. One can
clearly see that the interpolation does not suffer from the discontinuity. The
interpolation error has block-like structures, which are consistent with the
binary coding.

We now demonstrate the integration of internal variables in the TT format.
Once we obtain the integrated TT, we can readily evaluate it at any $\tau$ and
flavor pair. Figure 2(d) shows the self-energy evaluated on a dense grid for
the typical components. Note that $\Sigma^{(2)}$ and $\Pi^{(2)}$ are much
smaller than $\Sigma^{(1)}$ and $\Pi^{(1)}$ in amplitude, justifying the
validity of the one-shot calculation of $\Sigma^{(2)}$ and $\Pi^{(2)}$. The
$\tau$ dependence of $\Sigma^{(2)}$ and $\Pi^{(2)}$ does not exhibit any
signature of the discretization due to the exponential resolution. We further
tested the correctness of the evaluated self-energies by projecting them onto
the IR basis: We confirmed that the expansion coefficient decays as fast as
expected (not shown).

In this letter, we studied the low-rank structures of tensor trains (TTs)
representing the integrands in the multidimensional integrals of self-energy
Feynman diagrams using a prototype multiorbital electron-phonon model. We
analyzed second-order diagrams for the electron self-energy and phonon self-
energy. We discretized the continuous imaginary time variable via the quantics
representation and embedded discrete degrees of freedom, such as orbitals and
vibrational modes, in the same tensor train. We discovered that the integrands
possess low-rank structures, allowing for efficient data compression and
exponential convergence of the interpolation error (Fig. 1), regardless of
symmetry breaking by a superconducting transition/an external field.

We revealed that the discontinuity in the integrands for the self-energy does
not affect the compactness of the quantics representation [Fig. 3] and the
exponential convergence of the interpolation error. This robustness removes
the need for partitioning the imaginary-time domain aligned with the
discontinuity surfaces in conventional numerical integration approaches.

An interesting future direction is investigating the properties of the
electron-phonon model by incorporating electron-electron interactions, which
were ignored for simplicity in this study. Technically, various types of
optimization are immediately possible, such as using a sparse grid for the
external imaginary-time variable [23, 24, 25, 32] or automatic optimization of
index order in TT and tensor-network topology following the idea in Ref. 33.
Although the exponential convergence observed in the present study is
promising, the convergence properties at higher-order expansion orders remain
to be investigated. Extending the present approach to replacing QMC sampling
at higher expansion orders is interesting. Developing a hybrid approach of
QTCI and tailor-made algorithms for specific types of Feynman diagrams, e.g.,
Ref. 34, is another interesting future research topic.

###### Acknowledgements.

H.S. was supported by JSPS KAKENHI Grants No. 21H01041, No. 21H01003, and No.
23H03817 as well as JST FOREST Grant No. JPMJFR2232, Japan. H.I. and H.S. were
supported by JST PRESTO Grant No. JPMJPR2012, Japan. N.O and S.H were
supported by JSPS KAKENHI Grants No. 21K03459 and No. 23H01130. H.S. Thanks,
Mark Ritter, Jan von Delft, Markus Wallerberger, and Jason Kaye, for the
fruitful discussions. We used SparseIR.jl [35] for self-consistent
calculations on the sparse frequency and time grids.

## References

  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ (Kluwer Academic/Plenum Publishers, New York, 2000).
  * Hirsch and Fye [1986] J. E. Hirsch and R. M. Fye, Monte carlo method for magnetic impurities in metals, Physical review letters 56, 2521 (1986).
  * Van Houcke  _et al._ [2010] K. Van Houcke, E. Kozik, N. Prokof’ev, and B. Svistunov, Diagrammatic monte carlo, Physics Procedia 6, 95 (2010).
  * Prokof’ev and Svistunov [2007] N. Prokof’ev and B. Svistunov, Bold diagrammatic monte carlo technique: When the sign problem is welcome, Physical review letters 99, 250201 (2007).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time monte carlo methods for quantum impurity models, Reviews of Modern Physics 83, 349 (2011).
  * Shinaoka  _et al._ [2015] H. Shinaoka, Y. Nomura, S. Biermann, M. Troyer, and P. Werner, Negative sign problem in continuous-time quantum monte carlo: Optimal choice of single-particle basis for impurity problems, Physical Review B 92, 195126 (2015).
  * Pan and Meng [2024] G. Pan and Z. Y. Meng, The sign problem in quantum monte carlo simulations, in  _Encyclopedia of Condensed Matter Physics_ (Elsevier, 2024) p. 879â893.
  * Hann  _et al._ [2017] C. T. Hann, E. Huffman, and S. Chandrasekharan, Solution to the sign problem in a frustrated quantum impurity model, Annals of Physics 376, 63 (2017).
  * Zhang  _et al._ [2022] X. Zhang, G. Pan, X. Y. Xu, and Z. Y. Meng, Fermion sign bounds theory in quantum monte carlo simulation, Physical Review B 106, 035121 (2022).
  * Mondaini  _et al._ [2022] R. Mondaini, S. Tarat, and R. T. Scalettar, Quantum critical points and the sign problem, Science 375, 418 (2022).
  * Fernández  _et al._ [2022] Y. N. Fernández, M. Jeannin, P. T. Dumitrescu, T. Kloss, J. Kaye, O. Parcollet, and X. Waintal, Learning feynman diagrams with tensor trains, Physical Review X 12, 041018 (2022).
  * Erpenbeck  _et al._ [2023] A. Erpenbeck, W.-T. Lin, T. Blommel, L. Zhang, S. Iskakov, L. Bernheimer, Y. Núñez-Fernández, G. Cohen, O. Parcollet, X. Waintal, and E. Gull, Tensor train continuous time solver for quantum impurity models, Phys. Rev. B Condens. Matter 107, 245135 (2023).
  * Gourianov  _et al._ [2022] N. Gourianov, M. Lubasch, S. Dolgov, Q. Y. van den Berg, H. Babaee, P. Givi, M. Kiffner, and D. Jaksch, A quantum-inspired approach to exploit turbulence structures, Nature Computational Science 2, 30 (2022).
  * Gourianov [2022] N. Gourianov,  _Exploiting the structure of turbulence with tensor networks_ , Ph.D. thesis, University of Oxford (2022).
  * Ye and Loureiro [2022] E. Ye and N. F. G. Loureiro, Quantum-inspired method for solving the Vlasov-Poisson equations, Phys. Rev. E 106, 035208 (2022).
  * Shinaoka  _et al._ [2023] H. Shinaoka, M. Wallerberger, Y. Murakami, K. Nogaki, R. Sakurai, P. Werner, and A. Kauch, Multiscale space-time ansatz for correlation functions of quantum systems based on quantics tensor trains, Physical Review X 13, 021015 (2023).
  * Ritter  _et al._ [2024] M. K. Ritter, Y. N. Fernández, M. Wallerberger, J. von Delft, H. Shinaoka, and X. Waintal, Quantics tensor cross interpolation for high-resolution parsimonious representations of multivariate functions, Physical Review Letters 132, 056501 (2024).
  * Takabayashi and Prassides [2016] Y. Takabayashi and K. Prassides, Unconventional high-t c superconductivity in fullerides, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 20150320 (2016).
  * Nomura  _et al._ [2016] Y. Nomura, S. Sakai, M. Capone, and R. Arita, Exotic s-wave superconductivity in alkali-doped fullerides, Journal of Physics: Condensed Matter 28, 153001 (2016).
  * Capone  _et al._ [2009] M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti, Colloquium: Modeling the unconventional superconducting properties of expanded a 3 c 60 fullerides, Reviews of Modern Physics 81, 943 (2009).
  * Gunnarsson [1997] O. Gunnarsson, Superconductivity in fullerides, Reviews of modern physics 69, 575 (1997).
  * Kaga  _et al._ [2022] Y. Kaga, P. Werner, and S. Hoshino, Eliashberg theory of the jahn-teller-hubbard model, Physical Review B 105, 214516 (2022).
  * Shinaoka  _et al._ [2022] H. Shinaoka, N. Chikano, E. Gull, J. Li, T. Nomoto, J. Otsuki, M. Wallerberger, T. Wang, and K. Yoshimi, Efficient ab initio many-body calculations based on sparse modeling of matsubara green’s function, SciPost Phys. Lect. Notes , 063 (2022).
  * Li  _et al._ [2020] J. Li, M. Wallerberger, N. Chikano, C.-N. Yeh, E. Gull, and H. Shinaoka, Sparse sampling approach to efficient ab initio calculations at finite temperature, Phys. Rev. B 101, 035144 (2020).
  * Shinaoka  _et al._ [2017] H. Shinaoka, J. Otsuki, M. Ohzeki, and K. Yoshimi, Compressing green’s function using intermediate representation between imaginary-time and real-frequency domains, Physical Review B 96, 035147 (2017).
  * Oseledets [2009] I. Oseledets, Approximation of matrices with logarithmic number of parameters, in  _Doklady Mathematics_, Vol. 80 (Springer, 2009) pp. 653–654.
  * Khoromskij [2011] B. N. Khoromskij, O (d log n)-quantics approximation of n-d tensors in high-dimensional numerical modeling, Constructive Approximation 34, 257 (2011).
  * Takahashi  _et al._ [2024] H. Takahashi, R. Sakurai, and H. Shinaoka, Compactness of quantics tensor train representations of local imaginary-time propagators (2024), arXiv:2403.09161 [cond-mat.str-el] .
  * Núñez Fernández  _et al._ [2024] Y. Núñez Fernández, M. K. Ritter, M. Jeannin, J.-W. Li, T. Kloss, O. Parcollet, J. von Delft, H. Shinaoka, and X. Waintal, (2024), in preparation.
  * [30] See Supplemental Material at URL-will-be-inserted-by-publisher for the numerical details of our simulations.
  * Note [1] We must include the factor $\beta^{2}$ due to changing variables: $x^{\prime}=\tau^{\prime}/\beta$ and $x^{\prime\prime}=\tau^{\prime\prime}/\beta$.
  * Kaye  _et al._ [2022] J. Kaye, K. Chen, and O. Parcollet, Discrete lehmann representation of imaginary time green’s functions, Phys. Rev. B Condens. Matter 105, 235115 (2022).
  * Hikihara  _et al._ [2023] T. Hikihara, H. Ueda, K. Okunishi, K. Harada, and T. Nishino, Automatic structural optimization of tree tensor networks, Phys. Rev. Res. 5, 013031 (2023).
  * Kaye  _et al._ [2023] J. Kaye, H. U. R. Strand, and D. GoleÅ¾, Decomposing imaginary time feynman diagrams using separable basis functions: Anderson impurity model strong coupling expansion (2023), arXiv:2307.08566 [cond-mat.str-el] .
  * Wallerberger  _et al._ [2023] M. Wallerberger, S. Badr, S. Hoshino, S. Huber, F. Kakizawa, T. Koretsune, Y. Nagai, K. Nogaki, T. Nomoto, H. Mori, _et al._ , sparse-ir: Optimal compression and sparse sampling of many-body propagators, SoftwareX 21, 101266 (2023).


* Shukla et al. (2016) Jainendra Shukla, Miguel Barreda-Ángeles, Joan Oliver, and Domenec Puig.  Muderi: Multimodal database for emotion recognition among intellectually disabled individuals.  In _The Eight International Conference on Social Robotics_ , 2016\.
  * Silva et al. (2018) Michel Silva, Washington Ramos, Joao Ferreira, Felipe Chamone, Mario Campos, and Erickson R Nascimento.  A weighted sparse sampling and smoothing frame transition approach for semantic fast-forward first-person videos.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 2383–2392, 2018.
  * Soomro et al. (2012) Khurram Soomro, Amir Roshan Zamir, and Mubarak Shah.  Ucf101: A dataset of 101 human actions classes from videos in the wild.  _arXiv preprint arXiv:1212.0402_ , 2012.
  * Thornton et al. (2019) Mark A Thornton, Miriam E Weaverdyck, and Diana I Tamir.  The social brain automatically predicts others’ future mental states.  _Journal of Neuroscience_ , 39(1):140–148, 2019\.
  * Türker et al. (2018) Bekir Berker Türker, Engin Erzin, Yücel Yemez, and Metin Sezgin.  Audio-visual prediction of head-nod and turn-taking events in dyadic interactions.  _Interspeech_ , pages 1741–1745, 2018.
  * Ueno et al. (2020) Ryosuke Ueno, Yukiko I. Nakano, Jie Zeng, and Fumio Nihei.  Estimating the intensity of facial expressions accompanying feedback responses in multiparty video-mediated communication.  _ICMI 2020 - Proceedings of the 2020 International Conference on Multimodal Interaction_ , 20:144–152, 10 2020.
  * van Doorn (2018) Felix van Doorn.  Rituals of leaving: Predictive modelling of leaving behaviour in conversation, 2018.
  * van Son et al. (2008) Rob J.J.H. van Son, Wieneke Wesseling, Eric Sanders, and Henk van den Heuvel.  The ifadv corpus: a free dialog video corpus.  In _Proceedings of the Sixth International Conference on Language Resources and Evaluation (LREC’08)_ , 2008.
  * Vandeventer et al. (2015) Jason Vandeventer, Andrew J. Aubrey, Paul L. Rosin, and A. David Marshall.  4d cardiff conversation database (4d ccdb): a 4d database of natural, dyadic conversations.  _Auditory-Visual Speech Processing_ , 2015.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In _Advances in neural information processing systems_ , pages 5998–6008, 2017.
  * Vella and Paggio (2013) Alexandra Vella and Patrizia Paggio.  Overlaps in maltese: a comparison between map task dialogues and multimodal conversational data.  In _4th Nordic Symposium on Multimodal Communication_ , pages 21–29, 2013.
  * Vinciarelli et al. (2009) Alessandro Vinciarelli, Maja Pantic, and Hervé Bourlard.  Social signal processing: Survey of an emerging domain.  _Image and vision computing_ , 27(12):1743–1759, 2009.
  * von Marcard et al. (2018) Timo von Marcard, Roberto Henschel, Michael J. Black, Bodo Rosenhahn, and Gerard Pons-Moll.  Recovering accurate 3d human pose in the wild using imus and a moving camera.  _European Conference on Computer Vision (ECCV)_ , pages 601–617, 2018\.
  * Walker et al. (2017) Jacob Walker, Kenneth Marino, Abhinav Gupta, and Martial Hebert.  The pose knows: Video forecasting by generating pose futures.  In _Proceedings of the IEEE international conference on computer vision_ , pages 3332–3341, 2017.
  * Walsh et al. (2020) Kevin S Walsh, David P McGovern, Andy Clark, and Redmond G O’Connell.  Evaluating the neurophysiological evidence for predictive processing as a model of perception.  _Annals of the new York Academy of Sciences_ , 1464(1):242–268, 2020.
  * Wang et al. (2021a) Chenxi Wang, Yunfeng Wang, Zixuan Huang, and Zhiwen Chen.  Simple baseline for single human motion forecasting.  _ICCV SoMoF Workshop_ , 2021a.
  * Wang et al. (2021b) Jiashun Wang, Huazhe Xu, Medhini Narasimhan, Xiaolong Wang, and Uc San Diego.  Multi-person 3d motion prediction with multi-range transformers.  _Thirty-Fifth Conference on Neural Information Processing Systems_ , 2021b.
  * Wang et al. (2004) Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE transactions on image processing_ , 13(4):600–612, 2004.
  * Wolfert et al. (2021a) Pieter Wolfert, Jeffrey M Girard, Taras Kucherenko, and Tony Belpaeme.  To rate or not to rate: Investigating evaluation methods for generated co-speech gestures.  In _Proceedings of the 2021 International Conference on Multimodal Interaction_ , pages 494–502, 2021a.
  * Wolfert et al. (2021b) Pieter Wolfert, Nicole Robinson, and Tony Belpaeme.  A review of evaluation practices of gesture generation in embodied conversational agents.  _arXiv preprint arXiv:2101.03769_ , 2021b.
  * Woo et al. (2021) Jieyeon Woo, Catherine Pelachaud, and Catherine Achard.  Creating an interactive human/agent loop using multimodal recurrent neural networks.  _WACAI_ , 2021.
  * Yang et al. (2018) Dingdong Yang, Seunghoon Hong, Yunseok Jang, Tianchen Zhao, and Honglak Lee.  Diversity-sensitive conditional generative adversarial networks.  In _International Conference on Learning Representations_ , 2018.
  * Yang et al. (2021) Fangkai Yang, Yuan Gao, Ruiyang Ma, Sahba Zojaji, Ginevra Castellano, and Christopher E. Peters.  A dataset of human and robot approach behaviors into small free-standing conversational groups.  _PLoS ONE_ , 16, 2021.
  * Yasar and Iqbal (2021) Mohammad Samin Yasar and Tariq Iqbal.  A scalable approach to predict multi-agent motion for human-robot collaboration.  _IEEE Robotics and Automation Letters_ , 6(2):1686–1693, 2021.
  * Yonetani et al. (2016) Ryo Yonetani, Kris M. Kitani, and Yoichi Sato.  Recognizing micro-actions and reactions from paired egocentric videos.  _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 2629–2638, 2016.
  * Yuan and Kitani (2020) Ye Yuan and Kris Kitani.  Dlow: Diversifying latent flows for diverse human motion prediction.  In _European Conference on Computer Vision_ , pages 346–364. Springer, 2020.
  * Yuan and Kitani (2019) Ye Yuan and Kris M Kitani.  Diverse trajectory forecasting with determinantal point processes.  In _International Conference on Learning Representations_ , 2019.
  * Zhang and Yang (2021) Yu Zhang and Qiang Yang.  A survey on multi-task learning.  _IEEE Transactions on Knowledge and Data Engineering_ , 2021.
  * Zhao et al. (2019) Hang Zhao, Zhicheng Yan, Heng Wang, and Lorenzo Torresani.  Hacs: Human action clips and segments dataset for recognition and temporal localization.  _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ , pages 8667–8677, 2019.
  * Zhao and Dou (2020) Yang Zhao and Yong Dou.  Pose-forecasting aided human video prediction with graph convolutional networks.  _IEEE Access_ , 8:147256–147264, 2020.