# Sparse Box-fitting Least Squares

Aviad Panahi School of Physics and Astronomy, Raymond and Beverly Sackler
Faculty of Exact Sciences, Tel Aviv University, Tel Aviv 6997801, Israel  Shay
Zucker Porter School of the Environment and Earth Sciences, Raymond and
Beverly Sackler Faculty of Exact Sciences, Tel Aviv University, Tel Aviv
6997801, Israel

###### Abstract

We present a new implementation of the commonly used Box-fitting Least Squares
(BLS) algorithm, for the detection of transiting exoplanets in photometric
data. Unlike BLS, our new implementation – Sparse BLS (SBLS), does not use
binning of the data into phase bins, nor does it use any kind of phase grid.
Thus, its detection efficiency does not depend on the transit phase, and is
therefore slightly better than that of BLS. For sparse data, it is also
significantly faster than BLS. It is therefore perfectly suitable for large
photometric surveys producing unevenly-sampled sparse light curves, such as
Gaia.

Astronomy data analysis – Transit photometry – Exoplanet detection methods –
Algorithms

##  1 Introduction

The detection of transiting exoplanets is a very active research topic, with
implications ranging from physics to astrobiology. Since the first discovery
of a transiting exoplanet (Charbonneau et al., 2000), various approaches have
been applied to detect the telltale signals of exoplanetary transits in the
observed stellar light curves. These approaches mostly focus on the detection
of periodic dimmings in the light curves (e.g. Defaÿ et al., 2001; Kovács et
al., 2002; Renner et al., 2008; Hippke & Heller, 2019). The Box-fitting Least
Squares (BLS) algorithm (Kovács et al., 2002, hereafter KZM) is the most
widely-used approach for the detection of periodically transiting exoplanets.
It is based on a simplified box-shaped model of a strictly periodic transit,
which is characterized by only five parameters:

  * •

$P$ Period

  * •

$T_{0}$ Epoch of mid-transit

  * •

$w$ Transit duration

  * •

$M$ Out-of-transit mean magnitude

  * •

$d$ Transit depth

This simple approximation of a transit shape requires less parameters than a
more accurate one (e.g. Mandel & Agol, 2002; Hippke & Heller, 2019).
Furthermore, the two last parameters, $M$ and $d$, are determined analytically
by the data and the temporal parameters $P$, $T_{0}$ and $w$, and thus the
search space is three-dimensional and includes only the first three
parameters. Therefore, by the principle of parsimony (’Occam’s Razor’), it
brings about some statistical advantages (e.g. Jefferys & Berger, 1992),
especially when the data are scarce and their quality does not justify a more
detailed model. In addition to the statistical aspects, the simplicity also
translates to an algorithmic advantage, by requiring a shorter computation
time. In large surveys, with many light curves, using an overly complicated
model can result in prohibitively long run times. Thus, in the case of low-
cadence surveys such as Gaia (Gaia Collaboration et al., 2016), fitting
detailed models to each light curve, for the purpose of detection, is even
likely to be counterproductive, in both efficiency and run time. It is
therefore necessary to come up with an adaptation of BLS specifically designed
for such cases, of sparse light curves.

The popularity of BLS resulted in a variety of implementations and
improvements (e.g. Grziwa et al., 2012; Ofir, 2014; Hippke & Heller, 2019).
Many of these implementations are based on the original FORTRAN BLS code, in
which the light curve is phase-folded according to each trial period and
binned in phase. A double iteration is then performed over the bins to test
different beginning (ingress) and end (egress) phases of the transit.
Equivalently, some implementations (e.g. Collier Cameron et al., 2006) scan a
different parameter space, replacing the double iteration over phase bins with
iteration over discrete grids of phase and duration. Binning, or alternatively
discrete phase and duration grids, are useful when dealing with high-cadence
photometry with thousands of samples, mainly in terms of computing time.
However, they necessarily lose part of the information in the data, such as
the exact phases of all the samples. In the case of low-cadence light curves
containing only hundreds of samples, in which information is more scarce and
computing times are shorter, they might only introduce loss of information,
without any significant improvement in running time.

In this paper we present a new BLS implementation that does not rely on
binning nor on phase and duration grids. In Sec. 2 we present and detail this
novel implementation, and in Sec. 3 we discuss a few implementation issues
(such as run-time complexity). Section 4 discusses the performance of our
implementation, and we conclude in Sec. 5 with some final comments.

##  2 Sparse BLS

The implementation we introduce here, Sparse-BLS111A Java code is provided in:
https://github.com/aviadpi/SparseBLS (SBLS), essentially follows the same
original approach described in KZM, except the data is left unbinned. The
essence of the BLS approach is the identification of two brightness levels, in
and out of transit, and estimate them in order to compute a goodness-of-fit
statistic (e.g. $\chi^{2}$). When scanning trial configurations of $P$,
$T_{0}$ and $w$, only configuration changes that affect the inclusion of
points in the transit affect the goodness of fit, mainly by modifying the
estimates of the two levels. A naive grid scanning scheme that ignores the
phases of the light-curve samples might thus perform unnecessary computations
when the configuration change has no practical effect. Alternatively,
depending on the exact sample times, it might miss other configuration changes
that take place between grid points and do have an effect. SBLS, on the other
hand, focuses only on configuration changes that affect the inclusion of
points in the transit. In fact, for a given period, SBLS effectively scans the
space of transit phases and durations continuously, by checking only the
configurations in which the goodness-of-fit changes due to inclusion and
exclusion of points in transit.

SBLS scans a single grid of trial periods, followed by a double iteration over
indices in the phase-folded light curve. However, the indices now refer to the
_data points_ themselves, instead of the phase bins. This requires, of course,
sorting the points by phase for each trial period. Once the light curve is
sorted by phase, the double iteration looks for the points in transit. We
iterate over all data points for the first in-transit point, with index
$i_{1}$. For the last in-transit point, with index $i_{2}$, we only need to
scan points immediately following $i_{1}$, up to some predetermined maximum
duration (e.g. Ofir, 2014). We recommend setting a lower limit to the number
of points in transit, to reduce the chance of false detections caused by
spurious outliers. This of course depends on the total number of samples
available. For example, in our test runs (see Sec. 4), we simulated light
curves with a very small number of data points – $100$. We have set a lower
limit of three to the number of points in transits. This ensured a reasonable
performance while avoiding a substantial increase in spurious detections.

We repeat here the derivations presented in KZM. We denote the original time
series indices by $j$, while using $i$ for the indices of the time series
sorted by phase (’phase-folded indices’). Given a time series of values
$x_{j}$ and errors $\sigma_{j}$, we first convert it to a zero-weighted-mean
dataset with normalized weights $\left(\tilde{x_{j}},\tilde{w_{j}}\right)$:

$\displaystyle w_{j}$ $\displaystyle=\sigma_{j}^{-2}$ $\displaystyle W$
$\displaystyle=\sum\limits_{j=1}^{N}w_{j}$ $\displaystyle\tilde{w_{j}}$
$\displaystyle=\dfrac{w_{j}}{W}$ (1) $\displaystyle\mu$
$\displaystyle=\sum\limits_{j=1}^{N}\tilde{w_{j}}x_{j}$
$\displaystyle\tilde{x_{j}}$ $\displaystyle=x_{j}-\mu$ (2)

For each trial period $P$ we compute the phase values of the folded time
series using the modulo operation222We assume that the modulo operation
returns positive values for negative dividends. This is not necessarily the
case in all programming languages.:

$\phi_{i}=\frac{t_{i}\mod P}{P}$ (3)

For a given set of in-transit phase-folded indices
$\left\\{i_{1},i_{1}+1,\ldots,i_{2}\right\\}$, we calculate the Signal Residue
(SR):

$\displaystyle s$
$\displaystyle=\sum\limits_{i=i_{1}}^{i_{2}}\tilde{w_{i}}\tilde{x_{i}}$
$\displaystyle r$ $\displaystyle=\sum\limits_{i=i_{1}}^{i_{2}}\tilde{w_{i}}$
$\displaystyle\mathrm{\textit{SR}}=\left(\dfrac{s^{2}}{r(1-r)}\right)^{1/2}$
(4)

The $s$ and $r$ values can be used to estimate the stellar magnitude
$\left(M\right)$ and the transit depth $\left(d\right)$:

$\displaystyle\hat{M}$ $\displaystyle=\mu-\frac{s}{1-r}$
$\displaystyle\hat{d}$ $\displaystyle=\frac{s}{r\left(1-r\right)}$ (5)

The maximal SR value is saved for each trial period, resulting in a figure-of-
merit assigned to each period, i.e. a periodogram.

The calculation of the $s$ and $r$ values in SBLS makes use of memoization,
which reduces the computing time significantly. Memoization is a common
technique to optimize code execution by caching and retrieving interim values
of various variables. SBLS benefits from a memoization method we dubbed
”sliding window”, in which the values of $s$ and $r$ are sequentially updated
only once in each loop iteration, by subtracting or adding the necessary value
for every incremental change of either $i_{1}$ or $i_{2}$.

Figure 1 shows the BLS and SBLS spectra obtained for a simulated test case of
a low-cadence light curve with $200$ samples over a baseline period of
$1\,000$ days and a noise standard deviation of $1$ mmag. The simulated
transit parameters are based on those of HD$209458$b , with a period of
$3.52472$ days, transit duration of $0.128$ days and a depth of $0.0164$
magnitude (Charbonneau et al., 2000).

Figure 1 demonstrates another feature of SBLS: although the prominent peak in
the two spectra is very similar, there is a slight difference manifested as a
small slope in the BLS spectrum. This is related to the number of transit
configurations tested for each trial period, which is proportional to the
period in the specific phase grid scheme we used, and thus fewer
configurations are tested in the higher frequencies. The SBLS spectrum does
not use any scheme for binning or phase grid, and therefore the number of
configurations tested in each period depends only on the number of data
points.

Figure 1: Spectra of BLS (Top) and SBLS (Bottom). See text for details of the
simulated signal. In both spectra the most prominent peak is at the correct
frequency.

##  3 Implementation issues

###  3.1 Run-time complexity

Let us compare the estimated time complexity of the BLS and SBLS
implementations. We can safely omit the identical preparatory stages in which
we subtract the mean and normalize, especially since their complexity is only
$\mathcal{O}(N)$ operations. For the comparison we chose to use the equivalent
BLS implementation (e.g. Collier Cameron et al., 2006), rather than the
implementation that uses binning. This implementation scans a three-
dimensional parameter space of periods, phases and durations, with grid sizes
of $N_{P},N_{T_{0}},N_{w}$, respectively. For each configuration BLS then
proceeds to calculate SR, an $\mathcal{O}(N)$ procedure, resulting in a total
time complexity of $\mathcal{O}(N\cdot N_{P}\cdot N_{T_{0}}\cdot N_{w})$.
Thus, the BLS time complexity depends linearly on $N$ (assuming
$N_{P},N_{T_{0}},N_{w}$ are determined independently of $N$).

SBLS, on the other hand, scans a single grid of trial periods, and for each
period orders the samples by their phases, requiring $\mathcal{O}(N\log N)$
operations per trial period. Iteration over the in-transit indices
$\left(i_{1},i_{2}\right)$ requires $\mathcal{O}(N^{2})$ operations, though
this is only an upper limit (and therefore a worst-case estimate), as $i_{2}$
is usually not too far from $i_{1}$. The calculation of $SR$ would naively
require additional $\mathcal{O}(N)$ steps, but can be reduced to
$\mathcal{O}(1)$ using memoization. The total time complexity is therefore
$\mathcal{O}\left(N_{P}\cdot(N\log
N+N^{2})\right)\xrightarrow{}\mathcal{O}(N_{P}\cdot N^{2})$. In cases where
the phase and duration grid sizes, $N_{T_{0}}$ and $N_{w}$, depend on $N$, the
difference between BLS and SBLS time complexities would only increase, making
BLS less preferable over SBLS for large $N$. However, as we show in Sec. 4,
the balance might shift only at quite large values of $N$, and therefore SBLS
still has superior performance for small $N$.

###  3.2 Duration of transit

In this paper, phase values, denoted by $\phi$, are defined so that they are
dimensionless numbers ranging between $0$ and $1$. For any pair of transit-
scanning indices in the phase-folded light curve, $\left(i_{1},i_{2}\right)$,
we find their out-of-transit neighbours $\left(i_{1}^{-},i_{2}^{+}\right)$ and
use them to estimate ingress and egress phase values, while accounting for
wrap-around phases:

$\displaystyle\quad\quad\quad\quad\quad i_{1}^{-}=\left(i_{1}-1\right)\mod N$
(6) $\displaystyle\quad\quad\quad\quad\quad i_{2}^{+}=\left(i_{2}+1\right)\mod
N$ (7)
$\displaystyle\phi^{\text{ingress}}=\begin{cases}\frac{\phi[i_{1}]+\phi[i_{1}^{-}]}{2}&\text{if\quad}\phi[i_{1}^{-}]<\phi[i_{1}]\\\
\\\
\frac{\phi[i_{1}]+\phi[i_{1}^{-}]}{2}-\frac{1}{2}&\text{if\quad}\phi[i_{1}^{-}]>\phi[i_{1}]\end{cases}$
(8)
$\displaystyle\phi^{\text{egress}}=\begin{cases}\frac{\phi[i_{2}]+\phi[i_{2}^{+}]}{2}&\text{if\quad}\phi[i_{2}]<\phi[i_{2}^{+}]\\\
\\\
\frac{\phi[i_{2}]+\phi[i_{2}^{+}]}{2}+\frac{1}{2}&\text{if\quad}\phi[i_{2}]>\phi[i_{2}^{+}]\end{cases}$
(9)

Figure 2 serves to illustrate schematically the calculations detailed above.

Figure 2: Schematic illustration of the phase calculations in Eqs. 6–9. The
dashed magenta lines mark the phases of the first and last samples in the
transit ($i_{1}$ and $i_{2}$), while the dashed blue lines mark the phases of
the neighbouring samples out of transit. The red dashed-dotted lines mark the
estimated phases of the ingress and egress obtained by averaging the phases of
the corresponding in-transit and out-of-transit phases.

Finally we estimate the transit duration in units of time, using
$\phi^{\text{ingress}}$ and $\phi^{\text{egress}}$:

$\displaystyle w$
$\displaystyle=\left[\left(\phi^{\text{egress}}-\phi^{\text{ingress}}\right)\mod
1\right]\cdot P$ (10)

Figure 3 provides a pseudocode of the SBLS algorithm, depicting the steps
necessary for memoization, using the above definitions of the ingress and
egress phases, and the duration.

Input: Light curve LC$=\\{(t_{j},x_{j},\sigma_{j})\,|\,j=0\ldots N-1\\}$

Output: Periodogram, in-transit indices

Initialize LC;

for  _trial period $P$_ do

Sort LC by phase with $P$ (Eq. 3):
$\left(\phi_{i},x_{i}\right)\leftarrow\left(t_{j},x_{j}\right)$

${s},{r},{s}_{\text{prev}},{r}_{\text{prev}}\leftarrow 0$

for  _phase index $i_{1}\in[0,N-1]$ _ do

${w}\leftarrow{0}$

${s}\leftarrow{s}_{\text{prev}}$

${r}\leftarrow{r}_{\text{prev}}$

$i_{2}\leftarrow i_{2,\text{prev}}$

if  _$i_{1} >0$_ then

subtract values in $i_{1}^{-}$ and $i_{2}$ from ${s,r}$

estimate $\phi^{\text{ingress}}$ using Eq. 8

while  _$w <$ maxDuration_ do

if  _too few points in transit_ then

add values in $i_{2}$ to ${s,r}$

${s}_{\text{prev}}\leftarrow{s}$

${r}_{\text{prev}}\leftarrow{r}$

$i_{2,\text{prev}}\leftarrow i_{2}$

$i_{2}\leftarrow(i_{2}+1)\mod N$

continue

add values in $i_{2}$ to ${s,r}$

estimate $\phi^{\text{egress}}$ using Eq. 9

estimate $w$ using Eq. 10

calculate SR using ${s,r}$

if  _SR $>$ MaxSR_ then

MaxSR $\leftarrow$ SR

save in-transit indices

$i_{2}\leftarrow(i_{2}+1)\mod N$

save MaxSR to Periodogram

MaxSR $\leftarrow 0$

return _Periodogram , in-transit indices_

Figure 3: SBLS algorithm

###  3.3 Nearly sorted time series

In cases where the frequency grid is very dense, e.g. frequency steps of
$10^{-5}\,\text{~{}d}^{-1}$ with a time baseline of $T\sim 1\,000$ d, the
phase-folded order of the samples changes only slightly for successive trial
frequencies. In such cases, the sorting time complexity can be reduced by
applying sorting techniques tailored for nearly-sorted arrays on the light
curve which is already partially sorted from the previous trial frequency
(e.g. insertion sort; Cook & Kim, 1980). However, since the conventional
sorting time complexity, $\mathcal{O}(N\log N)$, is negligible compared to the
time spent in the double iteration, $\mathcal{O}(N^{2})$, the improvement in
run time is only on the order of a few percent of the total run time (less
than $5\%$ in our tests).

##  4 Performance

###  4.1 Detection efficiency

As KZM have already shown, binning has the immediate effect of introducing
some dependence of the performance on the transit phase, since the transit
phase need not necessarily align with the arbitrary binning scheme. This is
also the case for algorithms that scan discrete phase and duration grids. By
avoiding arbitrary binning or phase grids, the performance of SBLS is by
definition independent of the transit phase.

We present here one of the tests we have performed for comparing the
performance of BLS and SBLS, by applying them on two separate datasets – one
containing pure white Gaussian noise without any signal, and another
containing randomly generated transit signals at a predetermined signal-to-
noise ratio (SNR) as defined by KZM. For each dataset we simulated $1\,000$
light curves, each containing $100$ points with uniformly-distributed random
sampling times in the range of $\left[0,1\,000\right]$ days. We generated
transit signals with a transit depth of $d=0.01$ mag, noise with standard
deviation of $\sigma=3$ mmag, using frequencies drawn from a log-uniform
distribution in the range $\left[0.1,2.0\right]\mathrm{~{}d}^{-1}$, transit
phases drawn from a uniform distribution in the range $\left[0,1\right]$ and
transit duration set to $3\%$ of the drawn period (shorter duration might not
be detectable at all either by BLS or by SBLS with only $100$ samples). Both
BLS and SBLS scanned frequencies in the same range as the simulated
frequencies, with a frequency step of $10^{-4}\mathrm{~{}d}^{-1}$. For the BLS
we used a phase grid size of $N_{T_{0}}=1\,000$ and $N_{w}\leq 16$ for the
transit duration grid, ranging between $0.02\mathrm{~{}d}$ and the upper limit
for the transit duration, used in both implementations. The upper limit we
used for the duration was proportional to $P^{1/3}$, similarly to that used by
Collier Cameron et al. (2006). We used the maximum $SR$ value as the detection
statistic.

We chose to ignore the correctness of the frequency where the maximum $SR$ was
obtained, since if there were only few points in transit (due to the sparse
sampling), there could be several equivalently plausible frequencies that led
to the same affiliation of points to transit. In real life this issue should
be resolved by follow-up observations. Thus, a true detection was declared
when the returned maximum $SR$ value exceeded the threshold value when applied
on light curves containing a signal, disregarding its frequency.

We ended up with four distributions for the test statistic – ”Noise” and
”Signal” for each one of the two implementation. Figure 4 presents histograms
of the four distributions. Apparently, the two implementations produced very
similar distributions, as expected. In both of them there is considerable
overlap between the Signal and Noise histograms. This is related to the
relatively difficult detection challenge – only $100$ samples and considerable
noise. When the number of samples increases, the two distributions are
presenting diminishing overlap.

The Noise and Signal distributions described above can be used to compile a
Receiver Operating Characteristic (ROC) curve. A ROC curve (e.g. Fawcett,
2006) is a graphical tool to illustrate the performance of a detection scheme,
and is widely used in signal detection theory and related disciplines. As the
detection threshold increases, the rate of positive detection decreases – both
true (correctly identifying a transit signal as such) and false (wrongly
identifying a pure noise light curve as one that contains transit). The ROC
curve presents the way the true positive rate and the false positive rate are
related, by scanning different detection thresholds and counting the number of
true and false detections above them.

Figure 4: Sample distributions of the Maximum SR values for datasets
containing white noise and transit signals, of BLS (top) and Sparse-BLS
(bottom).

Figure 5 presents the two ROC curves resulting from the distributions in Fig.
4. The upper panel presents the two ROC curves, while the lower panel shows
the difference between them. It is clearly seen that the SBLS curve is
dominating the BLS curve, i.e. the SBLS true positive rate is larger than that
of the BLS for almost every assumed false positive rate. This is a
manifestation of the non-sensitivity of the SBLS performance to the transit
phase. Of course, using a much finer grid for the transit phase and duration
may eventually improve the performance of BLS, but that would come on the
expense of running time as was shown above in Sec. 3.1.

Figure 5: Top: Comparison of the ROC curves of BLS and SBLS. Bottom: The
difference between the two ROC curves. The details of the simulations used in
producing the curves are listed in the main text.

###  4.2 Run time comparison

We have implemented both BLS and SBLS in Java, and measured the run times of
both implementations for different sizes of datasets, using a $3.4$ GHz Intel
Octa-core i7 CPU. Figure 6 presents the results of this comparison. As
expected, one can see a linear dependence of the BLS run time on the size of
the light curve, and a quadratic dependence for SBLS. Indeed, for large light
curves SBLS is slower than BLS. However for sparse light curves, with less
than $5\,000$ points, the SBLS run time can be considerably shorter than that
of the BLS, even down to $2\%$.

Figure 6: Top: Run time of BLS and SBLS for light curves with up to $8\,000$
points. Bottom: Run time ratio: SBLS over BLS.

##  5 Conclusions

We have presented the SBLS (Sparse BLS) implementation of the BLS. Its main
advantage is that it is independent of any arbitrary phase grid or binning
scheme, unlike the original BLS of KZM. However, it is not a very efficient
implementation in terms of computing time, and therefore it might not be
suitable to long, high-cadence surveys, such as Kepler (Borucki et al., 2010)
or TESS (Ricker et al., 2014).

Nevertheless, the dependence of the complexity on the size of the light curve
is quadratic, and curiously enough, for data that is sparse enough, of less
than a few thousand samples per light curve, the linear dependence of the BLS
complexity is still considerably larger than the quadratic complexity for
SBLS. In those cases, the resource balance tilts towards SBLS, and it becomes
clear that SBLS is more efficient than BLS, both in terms of complexity and in
terms of detection efficiency. In fact, our complexity estimates of the BLS
runtime, which we presented in Sec. 3.1, were a little simplistic in assuming
that the number of bins (or alternatively grid cells) does not depend on the
light-curve size ($N$). If we introduce a dependence of $N_{T_{0}}$ and
$N_{w}$ on $N$, this would make the $N$-dependence of the BLS complexity
superlinear, increasing even more the breakeven point in the upper panel of
Fig. 6.

Thus, for large and sparse datasets such as Gaia, SBLS should be the tool of
choice for searching for periodic planetary transits.

This research was supported by the ISRAEL SCIENCE FOUNDATION (grant No.
848/16). We also acknowledge partial support by the Ministry of Science,
Technology and Space, Israel.

## References

  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Charbonneau et al. (2000) Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45, doi: 10.1086/312457
  * Collier Cameron et al. (2006) Collier Cameron, A., Pollacco, D., Street, R. A., et al. 2006, MNRAS, 373, 799, doi: 10.1111/j.1365-2966.2006.11074.x
  * Cook & Kim (1980) Cook, C. R., & Kim, D. J. 1980, Commun. ACM, 23, 620, doi: 10.1145/359024.359026
  * Defaÿ et al. (2001) Defaÿ, C., Deleuil, M., & Barge, P. 2001, A&A, 365, 330, doi: 10.1051/0004-6361:20000188
  * Fawcett (2006) Fawcett, T. 2006, Pattern Recognit. Lett., 27, 861 , doi: https://doi.org/10.1016/j.patrec.2005.10.010
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Grziwa et al. (2012) Grziwa, S., Pätzold, M., & Carone, L. 2012, MNRAS, 420, 1045, doi: 10.1111/j.1365-2966.2011.19970.x
  * Hippke & Heller (2019) Hippke, M., & Heller, R. 2019, A&A, 623, A39, doi: 10.1051/0004-6361/201834672
  * Jefferys & Berger (1992) Jefferys, W. H., & Berger, J. O. 1992, Am. Sci., 80, 64
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369, doi: 10.1051/0004-6361:20020802
  * Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, ApJ, 580, L171, doi: 10.1086/345520
  * Ofir (2014) Ofir, A. 2014, A&A, 561, A138, doi: 10.1051/0004-6361/201220860
  * Renner et al. (2008) Renner, S., Rauer, H., Erikson, A., et al. 2008, A&A, 492, 617, doi: 10.1051/0004-6361:200810148
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, 914320, doi: 10.1117/12.2063489


# Qubit-efficient quantum combinatorial optimization solver

Bhuvanesh Sundar<EMAIL_ADDRESS>Rigetti Computing, 775 Heinz Avenue,
Berkeley, California 94710, USA  Maxime Dupont  Rigetti Computing, 775 Heinz
Avenue, Berkeley, California 94710, USA

###### Abstract

Quantum optimization solvers typically rely on one-variable-to-one-qubit
mapping. However, the low qubit count on current quantum computers is a major
obstacle in competing against classical methods. Here, we develop a qubit-
efficient algorithm that overcomes this limitation by mapping a candidate bit
string solution to an entangled wave function of fewer qubits. We propose a
variational quantum circuit generalizing the quantum approximate optimization
ansatz (QAOA). Extremizing the ansatz for Sherrington-Kirkpatrick spin glass
problems, we show valuable properties such as the concentration of ansatz
parameters and derive performance guarantees. This approach could benefit
near-term intermediate-scale and future fault-tolerant small-scale quantum
devices.

An ongoing effort in quantum computing seeks to practically demonstrate
quantum utility for solving combinatorial optimization problems. Achieving
quantum utility could have a significant impact in areas which heavily utilize
optimization, such as logistics and supply chains. To this end, various
quantum algorithms have been developed and implemented, including annealing
approaches [1, 2], parameterized quantum algorithms with variational
optimization [3, 4, 5, 6], recursive techniques [7, 8, 9], and others inspired
by leading classical approaches [10, 11].

Quantum algorithms typically map each of the problem variables onto individual
qubits and recast the original optimization problem as the ground state search
for a Hamiltonian. One of the obstacles in demonstrating quantum utility is
the size of the problems one can address, given that current gate-based
quantum computers typically have tens to hundreds of qubits [12, 13, 14, 15,
16, 17, 18, 9, 19, 16], and quantum annealers have thousands of qubits [2]
111However, the effective problem sizes solved on quantum annealers are
smaller than the number of physical qubits due to embedding overheads
[choi2008minor].. Indeed, real-world problems typically have several thousands
of variables, which can be handled by classical solvers better than quantum
solvers today. A prevailing practical challenge in adding more qubits to the
computers is hardware noise, which limits their solution quality. In the
absence of quantum error correction [21], more qubits generally means larger
and therefore noisier quantum programs, compounding towards a deteriorating
output. Meanwhile, initial error-corrected quantum computers are expected to
have fewer logical qubits than necessary for quantum utility in optimization.
Therefore, qubit-efficient methods for optimization problems could benefit
near-term intermediate scale quantum devices as well as fault-tolerant small-
scale quantum computers [22, 23, 24, 25].

Several qubit-efficient methods have been proposed to overcome the size of
current quantum computers. Some of these are generic methods to break down a
large quantum circuit into an equivalent set of smaller quantum circuits that
can be executed on current devices [26, 27, 28, 29, 30], and have been used to
solve large optimization problems with more variables than qubits [31, 32]—at
the expense of large computational time overheads in classical or quantum
resources. Some approaches, specific to combinatorial optimization, employ
decompositions which stitch together an approximate solution for a large
optimization problem by solving smaller subproblems [33, 34, 35], and again
suffer from large overheads induced by the decomposition. Other strategies
employ nontrivial mappings between classical variables of the optimization
problem and qubits on the device, such as amplitude encoding [36, 37, 38],
correlation encoding [39], and multibasis encodings [40], each with some
challenges for quantum utility. For instance, the amplitude encoding strategy
is classically simulatable, while the other strategies use hardware-native
variational circuits which may be difficult to optimize [41, 42].

Here, we develop a qubit-efficient encoding based on a many-to-one mapping
from variables to qubits. In this approach, classical bitstrings are encoded
as entangled wave functions of fewer qubits. In addition, we physically
motivate a variational quantum circuit for this encoding, inspired by the
quantum approximate optimization algorithm (QAOA) [3]. While the QAOA uses
$q=N$ qubits to store $N$ variables [3], the proposed algorithm makes the
number of qubits $q$ a tunable parameter. Compared to the overhead of previous
decomposition methods [33, 34, 35, 26, 27, 28, 29, 30, 31, 32], this algorithm
only introduces a polynomial overhead in $Np$ over the QAOA for $p$ layers. We
show that this algorithm has several properties also present in the QAOA, such
as parameter clustering and performance guarantees for some classes of
problems [43, 44, 45, 46]. We numerically benchmark the algorithm for
Sherrington-Kirkpatrick (SK) spin glass models [47] and test it on Rigetti’s
AnkaaTM-9Q-3 superconducting chip.

Figure 1: (a) Qubit-efficient encoding of variables onto qubits via a many-to-
one mapping from variables to qubits. Graph nodes, which host variables, are
color-coded for convenience. All nodes with the same color are mapped to the
same data qubit (in grey). The value of the label qubits (purple) dictates
which node is stored in the data qubit. Links between nodes represent
interactions between nodes in the optimization problem [see text]. (b)
Parameterized quantum circuit to solve combinatorial optimization problems
with qubit-efficient encoding. The circuit consists of an initial layer of
Hadamard gates, followed by alternating layers of entangling
($e^{i\gamma_{p}\hat{\mathcal{H}}}$) and one-qubit gates ($U_{p}$). The
entangling gates in the $p$th layer is determined from measurements of the
wave function at the end of the $(p-1)$th layer. The parametric angles are
chosen classically to minimize the cost measured by the quantum computer.

_Qubit-efficient encoding.–_ Given $N$ binary variables
$\boldsymbol{z}=(z_{1},\cdots z_{N})$, we divide the variables into $N/d$
groups of $d$ variables each, as shown in Fig. 1(a), and denote each group
with a label $\ell$ which runs from 0 to $N/d-1$. Without loss of generality,
we assume $N$ is divisible by $d$. $d$ may be chosen by the user, for example
depending on the number of available qubits or the desired algorithm
performance as described later. The variables are stored in an entangled wave
function $\ket{\phi(\boldsymbol{z})}$ of $q=d+\log_{2}(N/d)$ qubits as

$\ket{\phi(\boldsymbol{z})}\equiv\sum\nolimits_{\ell=0}^{N/d-1}\lambda_{\ell}\ket{\ell}_{\rm
label}\otimes\ket{\boldsymbol{z}_{\ell}}_{\rm data},$ (1)

where $\ket{\ell}$ is the $\log_{2}(N/d)$-qubit product state in the binary
representation for $\ell$,
$\ket{\boldsymbol{z}_{\ell}}\equiv\ket{z_{(N/d-1)\ell+1},z_{(N/d-1)\ell+2},\cdots
z_{(N/d)\ell}}$ stores $d$ binary variables, and $\lambda_{\ell}$ are
arbitrary complex coefficients with $\sum_{\ell}|\lambda_{\ell}|^{2}=1$. This
wave function stores all $N$ binary variables of $\boldsymbol{z}$, and a
single projective measurement gives the values of $d$ variables,
$\boldsymbol{z}_{\ell}\equiv(z_{(N/d-1)\ell+1},z_{(N/d-1)\ell+2},\cdots
z_{(N/d)\ell})$, read from the data qubits conditioned on measuring the label
qubits in $\ell$. For convenience, we denote
$\ell_{i}\equiv\lfloor\frac{i}{d}\rfloor$ as the label value for the variable
$z_{i}$, and $d_{i}\equiv i\ \textrm{mod}\ d$ as the data qubit which stores
$z_{i}$ conditioned on the label qubits’ state being $\ket{\ell_{i}}$.

Figure 1(a) shows an example of the qubit-efficient encoding with $N=12,d=3$.
All variables $z_{i}$ with the same data qubit $d_{i}$ are given the same
color, for convenience of visualization. All variables within a blue ellipse
are given the same label $\ell$. In this case, the 12 variables are stored in
5 qubits.

In the case of $d=N$, $\ket{\phi(\boldsymbol{z})}\equiv\ket{z_{1},z_{2},\cdots
z_{N}}_{\textrm{d}ata}$ reduces to the usual way of storing binary variables
in a product state of the qubits. Nearly all the literature up to now for
solving optimization problems considered $d=N$ [4, 5]. For $d<N$, the wave
function stores the variables using a $(N/d)$-to-one map between the variables
and data qubits. This encoding describes a generic method to efficiently store
classical binary variables in qubits.

The case of $d=1$, first proposed in [48] and further explored in [49, 50, 51,
52] using a hardware-native ansatz, is notable since it requires only
$1+\log_{2}(N)$ qubits to store $N$ binary variables. However, we note that
quantum circuits in this case can be classically simulated in time scaling as
$O(Np)$ where $p$ parametrizes the circuit depth (see below), and similarly
for $d\lesssim\log_{2}N$. Under the widely held belief that quantum algorithms
may not be able to achieve exponential speedup for NP-hard problems, it is
unclear where the difficulty for the quantum algorithms at
$d\lesssim\log_{2}N$ lies, with potential candidates being a large circuit
depth $p$ or expensive computational time for parameter search.

The main aim of this paper is to develop a quantum algorithm that uses qubit-
efficient encoding to solve combinatorial optimization problems, and
investigate the cases with $1\leq d\leq N$.

_Binary optimization problem.–_ Binary optimization problems seek to minimize
a cost function $C(\boldsymbol{z})$ of variables $\boldsymbol{z}=(z_{1},\cdots
z_{N})$ with $z_{i}\in\\{\pm 1\\}$, where the function $C$ is generically
written as a polynomial in the variables, e.g.,
$C=\sum_{i}w_{i}z_{i}+\sum_{ij}w_{ij}z_{i}z_{j}+\cdots$. To develop a quantum
algorithm that minimizes this cost, one should express the cost as expectation
value(s) with respect to samples from the wave function $\ket{\psi}$. This
task is trivial for $d=N$, but nontrivial for $d<N$.

For simplicity, let us consider one of the singleton terms in the cost,
$w_{i}z_{i}$. Since the value of the variable $z_{i}$ can be obtained only
when the label bits are measured as $\ell_{i}$, an estimation
$\overline{z}_{i}$ for this variable is obtained from averaging the shot
values of the Pauli operator $\hat{Z}_{d_{i}}$, post-selecting only those
measurements with the label bits as $\ell_{i}$. This can be written as the
conditional expectation value,

$\bar{z}_{i}=\frac{\braket{\psi}{\hat{P}_{\ell_{i}}\hat{Z}_{d_{i}}}{\psi}}{\braket{\psi}{\hat{P}_{\ell_{i}}}{\psi}},$
(2)

where $\hat{P}_{\ell}=\ket{\ell}_{\textrm{label}}\bra{\ell}_{\textrm{label}}$
is a projection operator. Similarly, for variables $z_{i}$ and $z_{j}$ with
the same label, $\ell_{i}=\ell_{j}$, $z_{i}z_{j}$ can be estimated by post-
selecting on shots where the label bits were $\ell_{i}$, therefore the
estimate $\overline{z_{i}z_{j}}$ is

$\left(\overline{z_{i}z_{j}}\right)_{\ell_{i}=\ell_{j}}=\frac{\braket{\psi}{\hat{P}_{\ell_{i}}\hat{Z}_{d_{i}}\hat{Z}_{d_{j}}}{\psi}}{\braket{\psi}{\hat{P}_{\ell_{i}}}{\psi}}.$
(3)

If $\ell_{i}\neq\ell_{j}$, $z_{i}$ and $z_{j}$ are never measured together in
the same shot, therefore we estimate $\overline{z_{i}z_{j}}$ as
$\bar{z}_{i}\times\bar{z}_{j}$. In this way, quadratic terms of the form
$C(\boldsymbol{z})=\sum_{i<j}w_{ij}z_{i}z_{j}$, for example, can be expressed
as

$C[\psi]=\sum_{\begin{subarray}{c}i<j\\\
\ell_{i}=\ell_{j}\end{subarray}}w_{ij}\overline{z_{i}z_{j}}+\sum_{\begin{subarray}{c}i<j\\\
\ell_{i}\neq\ell_{j}\end{subarray}}w_{ij}\bar{z}_{i}\bar{z}_{j}.$ (4)

This construction can be straightforwardly extended to measure any polynomial
in $\boldsymbol{z}$.

The cost in Eq. (4) can be written as
$C[\psi]=\braket{\psi}{\hat{\mathcal{H}}[\psi]}{\psi}$ with

$\displaystyle\hat{\mathcal{H}}[\psi]=$
$\displaystyle\sum_{\begin{subarray}{c}i<j\\\
\ell_{i}=\ell_{j}\end{subarray}}w_{ij}\frac{\hat{P}_{\ell_{i}}\hat{Z}_{d_{i}}\hat{Z}_{d_{j}}}{\braket{\psi}{\hat{P}_{\ell_{i}}}{\psi}}$
$\displaystyle+$ $\displaystyle\frac{1}{2}\sum_{\begin{subarray}{c}i<j\\\
\ell_{i}\neq\ell_{j}\end{subarray}}w_{ij}\left(\bar{z}_{i}\frac{\hat{P}_{\ell_{j}}\hat{Z}_{d_{j}}}{\braket{\psi}{\hat{P}_{\ell_{j}}}{\psi}}+\bar{z}_{j}\frac{\hat{P}_{\ell_{i}}\hat{Z}_{d_{i}}}{\braket{\psi}{\hat{P}_{\ell_{i}}}{\psi}}\right).$
(5)

We highlight that the Hamiltonian $\hat{\mathcal{H}}[\psi]$ depends on
$\ket{\psi}$, a consequence of the cost being nonlinear in
$\ket{\psi}\bra{\psi}$. The first line in Eq. (Qubit-efficient quantum
combinatorial optimization solver) describes interactions between nodes with
the same group, and the second line in Eq. (Qubit-efficient quantum
combinatorial optimization solver) describes interactions between nodes in
different groups.

This construction of $\hat{\mathcal{H}}$ and $C$ smoothly connects in the
$d=N$ limit with their standard forms in the QAOA,
$\hat{\mathcal{H}}=\sum_{i}w_{i}\hat{Z}_{i}+\sum_{ij}w_{ij}\hat{Z}_{i}\hat{Z}_{j}+\cdots$
and $C=\braket{\hat{\mathcal{H}}}$. Therefore many known properties of the
QAOA may be expected to hold for values of $d$ close to $N$ as well, as
exemplified later.

Here, we restrict our attention to SK models [47], which have quadratic cost
functions of the form of Eq. (4) with fully connected random weights $w_{ij}$
precisely defined later. Optimization problems expressed as fully connected
graphs are ubiquitous in statistical mechanics and computer science [53].

_Ansatz.–_ The QAOA ansatz for extremizing $C$ is
$\ket{\psi_{p}}=e^{i\beta_{p}\hat{\mathcal{H}}_{x}}e^{i\gamma_{p}\hat{\mathcal{H}}}\ket{\psi_{p-1}}$
[3]. A natural extension of this ansatz generalizes $\hat{\mathcal{H}}$ to
$\hat{\mathcal{H}}[\psi_{p-1}]$. Here,
$\hat{\mathcal{H}}_{x}=\sum_{i=1}^{q}\hat{X}_{i}$ is the transverse field
mixer Hamiltonian, and $\hat{\mathcal{H}}[\psi_{p-1}]$ is the layer-dependent
cost Hamiltonian, which depends on the wavefunction $\ket{\psi_{p-1}}$ in the
previous layer.

The ansatz is drawn in Fig. 1(b). Determining the cost Hamiltonian in the
$p$th layer, $\hat{\mathcal{H}}[\psi_{p-1}]$, requires to first know $\bar{z}$
and $\braket{\psi_{p-1}}{\hat{P}_{\ell}}{\psi_{p-1}}$, which can be obtained
from projective measurements after $p-1$ layers. Circuit implementation of
$\exp(i\gamma\hat{\mathcal{H}})$ is discussed in [54]. Below, we exemplify
this ansatz using a SK problem instance, and argue for a slight modification
to the ansatz.

_N=4 example.–_

Figure 2: (a) Schematic of quantum circuit that implements a $p=1$ layer of
the ansatz, for a $N=4$ SK model instance at $d=2$. Orange and black dots are
$R_{x}$ and $R_{z}$ gates, and vertical lines are iSWAPs [54]. (b-c) Cost
landscape for the ansatz versus $\beta_{1}$ and $\gamma_{1}$, measured from
experiments on Rigetti’s AnkaaTM-9Q-3 chip, and numerically computed from an
ideal simulation of the circuit, with 10000 shots.

We consider a SK model problem with $N=4$ variables, where the weights
$w_{ij}$ are chosen randomly from $\\{\pm 1\\}$. The model weights for the
specific instance are given in [54]. For this example, we encode the problem
on 3 qubits using $d=2$. Figure 2(a) shows the circuit required to implement
the ansatz at $p=1$, expressed in terms of quantum gates native to Rigetti’s
AnkaaTM-9Q-3 superconducting quantum chip (see [54] for definitions of these
gates, and a circuit compilation of the ansatz).

Figures 2(b-c) plot the cost yielded by the ansatz versus $\beta_{1}$ and
$\gamma_{1}$, with results measured on Rigetti’s AnkaaTM-9Q-3 in (b), and
numerically obtained from a noiseless simulator in (c). We find extremely good
agreement between the landscapes, with a Pearson correlation coefficient of
0.997. There is a small $\sim 15\%$ scale reduction of the cost on
AnkaaTM-9Q-3, that is attributable to hardware noise. However, the $p=1$
ansatz does not produce the ground state cost (see [54]), and cannot produce
the ground state cost at any $p$ due to a symmetry explained below.

The SK problem instance has two degenerate solutions, $\boldsymbol{z}^{*}$ and
$-\boldsymbol{z}^{*}$, due to a global $\mathbb{Z}_{2}$ spin-flip symmetry.
Similarly, the ansatz also has a global spin-flip symmetry for the data
qubits. The ansatz at $p=1$ and the optimal angles [54] prepares an equal
superposition of the two ground state solutions $\boldsymbol{z}^{*}$ and
$-\boldsymbol{z}^{*}$, i.e.
$\ket{\psi_{p=1}}=(\ket{\phi(\boldsymbol{z}^{*})}+\ket{\phi(-\boldsymbol{z}^{*})})/\sqrt{2}$.
Yet, since $\bar{z}_{i}=0\ \forall i$ due to the spin-flip symmetry (see
[54]), the ansatz fails to produce the minimum cost. We note that this is
unique to the qubit-efficient encoding, and does not occur at $d=N$.

Thus, breaking the $\mathbb{Z}_{2}$ symmetry of the ansatz is needed to
produce the ground state cost. Therefore, we modify the ansatz to

$\ket{\psi_{p}}=e^{i\beta_{p}\hat{\mathcal{H}}_{x}}e^{i\gamma^{\prime}_{p}\hat{\mathcal{H}}_{z}}e^{i\gamma_{p}\hat{\mathcal{H}}[\psi_{p-1}]}\ket{\psi_{p-1}},$
(6)

where $\hat{\mathcal{H}}_{z}=\sum_{i\in\textrm{data}}\hat{Z}_{i}$ describes
one-body bias terms.

_Variational optimization.–_

Figure 3: Average approximation ratios produced by the optimal circuits in
Fig. 1(b), for ensemble of five Sherrington-Kirkpatrick spin glass model
problems with $N=64$ variables and $w_{ij}\in\\{\pm 1\\}$, at various values
of $d$.

We illustrate the performance of the symmetry-broken ansatz using an ensemble
of SK problem instances at $N=64$ and $128$ with $w_{ij}\in\\{\pm 1\\}$. We
optimize the $(\beta_{i},\gamma_{i},\gamma^{\prime}_{i})$ parameters in the
ansatz with a classical optimizer which seeks to minimize $C[\psi_{p}]$. We
use a basin-hopping optimizer [55] with L-BFGS-B as the local optimization
method [56]. To benchmark the cost produced by variational optimization, we
compare it against the ground state cost $C^{*}$ found with a classical
heuristic algorithm based on tabu search [57] implemented in MQLib [58].

In Fig. 3, we plot the optimized approximation ratios, $r=C/C^{*}$, averaged
over an ensemble of five SK problems with $N=64$ variables, using different
values of $d<N$. The dashed line at $d=N$ is plotted using values predicted in
Ref. [43] for large $N$. Our ansatz at $d<N,p=3$ produces results comparable
to the ansatz at $d=N,p=1$, despite having up to 9 times fewer qubits at
$d=1$. Similarly, the ansatz at $d=1,p=4$ produces results similar to the
ansatz at $d=N,p=2$.

In addition to the quality of the solution produced by an algorithm, the time
taken by the algorithm is also an important factor. A significant contributor
to the total time of variational quantum algorithms is the time taken to find
the optimal parameters. Figure 3 shows no clear distinction in the quality of
solution produced by different $d$ at the same $p$. But we argue next that
larger values of $d$ may require lesser parameter search time due to
clustering of optimal parameters for SK problem instances.

_Concentration of parameters.–_

Figure 4: Average approximation ratios for optimized parameters for five SK
problem instances (teal), and for 100 SK problem instances where the same
parameters were used for all the instances (yellow and magenta). Instances in
yellow have $w_{ij}=\pm 1$ while magenta instances have
$w_{ij}\sim\mathcal{N}(0,1)$. The parameters for the yellow and magenta curves
in (b) are related to those of (a) by a trivial rescaling factor [54].

Previous works have shown empirical and analytical evidence for clustering of
optimal parameters in the QAOA ($d=N$) for some classes of problems including
SK problems [59, 43, 60, 44, 61, 62]. It is interesting to ask whether a
similar property holds for the qubit-efficient ansatz. Figure 4 demonstrates
empirical evidence for parameter concentration for SK problems in the qubit-
efficient ansatz. We implemented the ansatz for solving 100 SK problem
instances, using the optimized parameters for one of the SK problem instances
from Fig. 3. The yellow dashed line in Fig. 4(a) plots the mean $r$ vs $p$ for
100 SK problem instances with $w_{ij}=\pm 1$, and shows that the mean $r$ for
100 problem instances closely follows the optimized $r$ for five problem
instances (solid teal curve). Furthermore, the magenta dashed line in Fig.
4(a), which plots the mean approximation ratio for problems with normally
distributed weights, $w_{ij}\sim\mathcal{N}(0,1)$, shows evidence that the
same parameters obtained from the case with $w_{ij}\in\\{\pm 1\\}$ can also be
used to solve the ensemble with $w_{ij}\sim\mathcal{N}(0,1)$. Our ansatz beats
a baseline performance that comes from a naive problem decomposition method,
explained in [54].

Figure 4(b) demonstrates another useful property that near-optimal parameters
for $N>64$ can be derived from those at $N=64$. We solve SK problems with
$N=128$, by reusing the same parameters from the $N=64$ case in (a) with a
trivial rescaling [54], and still find approximation ratios close to the
optimized $r$ for five problem instances (solid teal curve).

_Outlook.–_ We have shown that a QAOA-inspired algorithm for solving SK spin
glass problems, with a qubit-efficient encoding scheme which maps many
variables to each qubit in a register, achieves an approximation ratio within
$O(1)$ of that of the QAOA, while using substantially fewer qubits. We
demonstrated a proof-of-principle experimental benchmark of the algorithm for
solving a small SK model problem on Rigetti’s AnkaaTM-9Q-3 superconducting
chip. We derived a performance guarantee for our ansatz, and showed evidence
of parameter clustering at $d=4$, which could significantly alleviate
classical resources for parameter optimization. We did not find parameter
concentration at $d<4$, which could indicate one of the difficulties in
implementing the algorithm with $O(\log_{2}N)$ qubits, and therefore a
preference for larger $d$. However, this preference for more qubits should be
weighed against the increased hardware noise in the circuit when more qubits
are present. Moreover, one should consider the fact that the solution bit
string would be stored as an entangled final wave function which may be less
robust to noise when more qubits are used. Thus, investigating the optimal
value of $d$ on real hardware would be an interesting question to explore in
the future.

Our algorithm can be straightforwardly extended to other classes of
combinatorial problems, and is expected to have a trainable QAOA-like ansatz
that may have parameter clustering or at least good heuristics for seed
parameters, and lesser overheads than standard problem decomposition methods
[33, 34, 35, 26, 27, 28, 29, 30, 31, 32]. The qubit-efficient scheme may also
be applied to other use cases such as quantum machine learning, and is likely
to be useful in the near future as we enter the era of small-scale fault-
tolerant quantum computers [22, 23, 24, 25].

###### Acknowledgements.

We thank Stuart Hadfield, Mark Hodson, Marco Paini, Alex Place, Matt Reagor,
Graham Reid, and Davide Venturelli for valuable discussions. This work is
supported by the Defense Advanced Research Projects Agency (DARPA) under
Agreement No. HR-0011-23-3-0015. This research used resources of the National
Energy Research Scientific Computing Center, a DOE Office of Science User
Facility supported by the Office of Science of the U.S. Department of Energy
under Contract No. DE-AC02-05CH11231 using NERSC awards DDR-ERCAP0024427 and
ASCR-ERCAP0028951. The views, opinions and/or findings expressed are those of
the authors and should not be interpreted as representing the official views
or policies of the Department of Defense or the U.S. Government.

## References

  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum computation by adiabatic evolution, arXiv , 0001106 (2000).
  * King  _et al._ [2023] A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. Berkley, K. Boothby, S. Ejtemaee, C. Enderud, _et al._ , Quantum critical dynamics in a 5,000-qubit programmable spin glass, Nature 617, 61 (2023).
  * [3] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv , 1411.4028.
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nat. Rev. Phys. 3, 625 (2021).
  * Blekos  _et al._ [2024] K. Blekos, D. Brand, A. Ceschini, C.-H. Chou, R.-H. Li, K. Pandya, and A. Summer, A review on quantum approximate optimization algorithm and its variants, Phys. Rep. 1068, 1 (2024).
  * Zhu  _et al._ [2022] L. Zhu, H. L. Tang, G. S. Barron, F. A. Calderon-Vargas, N. J. Mayhall, E. Barnes, and S. E. Economou, Adaptive quantum approximate optimization algorithm for solving combinatorial problems on a quantum computer, Phys. Rev. Res. 4, 033029 (2022).
  * Bravyi  _et al._ [2020] S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, Obstacles to variational quantum optimization from symmetry protection, Phys. Rev. Lett. 125, 260505 (2020).
  * [8] L. T. Brady and S. Hadfield, Iterative quantum algorithms for maximum independent set: a tale of low-depth quantum algorithms, arXiv , 2309.13110.
  * Dupont  _et al._ [2023] M. Dupont, B. Evert, M. J. Hodson, B. Sundar, S. Jeffrey, Y. Yamaguchi, D. Feng, F. B. Maciejewski, S. Hadfield, M. S. Alam, _et al._ , Quantum-enhanced greedy combinatorial optimization solver, Sci. Adv. 9, eadi0487 (2023).
  * Dupont and Sundar [2024] M. Dupont and B. Sundar, Extending relax-and-round combinatorial optimization solvers with quantum correlations, Phys. Rev. A 109, 012429 (2024).
  * Egger  _et al._ [2021] D. J. Egger, J. Mareček, and S. Woerner, Warm-starting quantum optimization, Quantum 5, 479 (2021).
  * Harrigan  _et al._ [2021] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, _et al._ , Quantum approximate optimization of non-planar graph problems on a planar superconducting processor, Nat. Phys. 17, 332 (2021).
  * Ebadi  _et al._ [2022] S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J.-G. Liu, R. Samajdar, _et al._ , Quantum optimization of maximum independent set using rydberg atom arrays, Science 376, 1209 (2022).
  * [14] F. B. Maciejewski, S. Hadfield, B. Hall, M. Hodson, M. Dupont, B. Evert, J. Sud, M. S. Alam, Z. Wang, S. Jeffrey, _et al._ , Design and execution of quantum circuits using tens of superconducting qubits and thousands of gates for dense ising optimization problems, arXiv , 2308.12423.
  * Shaydulin  _et al._ [2024] R. Shaydulin, C. Li, S. Chakrabarti, M. DeCross, D. Herman, N. Kumar, J. Larson, D. Lykov, P. Minssen, Y. Sun, _et al._ , Evidence of scaling advantage for the quantum approximate optimization algorithm on a classically intractable problem, Sci. Adv. 10, eadm6761 (2024).
  * [16] N. Sachdeva, G. S. Harnett, S. Maity, S. Marsh, Y. Wang, A. Winick, R. Dougherty, D. Canuto, Y. Q. Chong, M. Hush, _et al._ , Quantum optimization using a 127-qubit gate-model ibm quantum computer can outperform quantum annealers for nontrivial binary optimization problems, arXiv , 2406.01743.
  * Moses  _et al._ [2023] S. A. Moses, C. H. Baldwin, M. S. Allman, R. Ancona, L. Ascarrunz, C. Barnes, J. Bartolotta, B. Bjork, P. Blanchard, M. Bohn, _et al._ , A race-track trapped-ion quantum processor, Phys. Rev. X 13, 041052 (2023).
  * Pelofske  _et al._ [2024] E. Pelofske, A. Bärtschi, and S. Eidenbenz, Short-depth qaoa circuits and quantum annealing on higher-order ising models, npj Quantum Inf. 10, 30 (2024).
  * [19] J. A. Montanez-Barrera and K. Michielsen, Towards a universal qaoa protocol: Evidence of quantum advantage in solving combinatorial optimization problems, arXiv , 2405.09169.
  * Note [1] However, the effective problem sizes solved on quantum annealers are smaller than the number of physical qubits due to embedding overheads [choi2008minor].
  * Roffe [2019] J. Roffe, Quantum error correction: an introductory guide, Contemp. Phys. 60, 226 (2019).
  * Bluvstein  _et al._ [2024] D. Bluvstein, S. J. Evered, A. A. Geim, S. H. Li, H. Zhou, T. Manovitz, S. Ebadi, M. Cain, M. Kalinowski, D. Hangleiter, _et al._ , Logical quantum processor based on reconfigurable atom arrays, Nature 626, 58 (2024).
  * AI [2023] G. Q. AI, Suppressing quantum errors by scaling a surface code logical qubit, Nature 614, 676 (2023).
  * [24] M. P. Da Silva, C. Ryan-Anderson, J. M. Bello-Rivas, A. Chernoguzov, J. M. Dreiling, C. Foltz, J. P. Gaebler, T. M. Gatterman, D. Hayes, N. Hewitt, _et al._ , Demonstration of logical qubits and repeated error correction with better-than-physical error rates, arXiv , 2404.02280.
  * Krinner  _et al._ [2022] S. Krinner, N. Lacroix, A. Remm, A. Di Paolo, E. Genois, C. Leroux, C. Hellings, S. Lazar, F. Swiadek, J. Herrmann, _et al._ , Realizing repeated quantum error correction in a distance-three surface code, Nature 605, 669 (2022).
  * Bravyi  _et al._ [2016] S. Bravyi, G. Smith, and J. A. Smolin, Trading classical and quantum computational resources, Phys. Rev. X 6, 021043 (2016).
  * Peng  _et al._ [2020] T. Peng, A. W. Harrow, M. Ozols, and X. Wu, Simulating large quantum circuits on a small quantum computer, Phys. Rev. Lett. 125, 150504 (2020).
  * Tang  _et al._ [2021] W. Tang, T. Tomesh, M. Suchara, J. Larson, and M. Martonosi, Cutqc: using small quantum computers for large quantum circuit evaluations, in  _Proceedings of the 26th ACM International conference on architectural support for programming languages and operating systems (New York)_ (2021) pp. 473–486.
  * DeCross  _et al._ [2023] M. DeCross, E. Chertkov, M. Kohagen, and M. Foss-Feig, Qubit-reuse compilation with mid-circuit measurement and reset, Phys. Rev. X 13, 041057 (2023).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, _et al._ , Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Bechtold  _et al._ [2023] M. Bechtold, J. Barzen, F. Leymann, A. Mandl, J. Obst, F. Truger, and B. Weder, Investigating the effect of circuit cutting in qaoa for the maxcut problem on nisq devices, Quantum Sci. Technol. 8, 045022 (2023).
  * [32] M. Dupont, B. Sundar, B. Evert, D. E. B. Neira, Z. Peng, S. Jeffrey, and M. J. Hodson, Quantum optimization for the maximum cut problem on a superconducting quantum computer, arXiv , 2404.17579.
  * Angone  _et al._ [2023] A. Angone, X. Liu, R. Shaydulin, and I. Safro, Hybrid quantum-classical multilevel approach for maximum cuts on graphs, in  _2023 IEEE High Performance Extreme Computing Conference (HPEC)_ (IEEE, 2023) pp. 1–7.
  * [34] B. Bach, J. Falla, and I. Safro, Mlqaoa: Graph learning accelerated hybrid quantum-classical multilevel qaoa, arXiv , 2404.14399.
  * [35] C.-Y. Liu and H.-S. Goan, Hybrid gate-based and annealing quantum computing for large-size ising problems, arXiv , 2208.03283.
  * Rančić [2023] M. J. Rančić, Noisy intermediate-scale quantum computing algorithm for solving an n-vertex maxcut problem with log (n) qubits, Phys. Rev. Res. 5, L012021 (2023).
  * [37] Y. Chatterjee, E. Bourreau, and M. J. Rančić, Solving various np-hard problems using exponentially fewer qubits on a quantum computer, arXiv , 2301.06978.
  * [38] Y. Tene-Cohen, T. Kelman, O. Lev, and A. Makmal, A variational qubit-efficient maxcut heuristic algorithm, arXiv , 2308.10383.
  * [39] M. Sciorilli, L. Borges, T. L. Patti, D. García-Martín, G. Camilo, A. Anandkumar, and L. Aolita, Towards large-scale quantum optimization solvers with few qubits, arXiv , 2401.09421.
  * Patti  _et al._ [2022] T. L. Patti, J. Kossaifi, A. Anandkumar, and S. F. Yelin, Variational quantum optimization with multibasis encodings, Phys. Rev. Res. 4, 033142 (2022).
  * Bittel and Kliesch [2021] L. Bittel and M. Kliesch, Training variational quantum algorithms is np-hard, Phys. Rev. Lett. 127, 120502 (2021).
  * [42] M. Larocca, S. Thanasilp, S. Wang, K. Sharma, J. Biamonte, P. J. Coles, L. Cincio, J. R. McClean, Z. Holmes, and M. Cerezo, A review of barren plateaus in variational quantum computing, arXiv , 2405.00781.
  * Farhi  _et al._ [2022] E. Farhi, J. Goldstone, S. Gutmann, and L. Zhou, The quantum approximate optimization algorithm and the sherrington-kirkpatrick model at infinite size, Quantum 6, 759 (2022).
  * [44] J. Wurtz and D. Lykov, The fixed angle conjecture for qaoa on regular maxcut graphs, arXiv , 2107.00677.
  * Wurtz and Love [2021] J. Wurtz and P. Love, Maxcut quantum approximate optimization algorithm performance guarantees for p¿ 1, Phys. Rev. A 103, 042612 (2021).
  * Basso  _et al._ [2022] J. Basso, E. Farhi, K. Marwaha, B. Villalonga, and L. Zhou, The quantum approximate optimization algorithm at high depth for maxcut on large-girth regular graphs and the sherrington-kirkpatrick model, in  _17th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2022)_ (2022) pp. 7:1–7:21.
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Solvable model of a spin-glass, Phys. Rev. Lett. 35, 1792 (1975).
  * Tan  _et al._ [2021] B. Y. L. Tan, M.-A. Lemonde, S. Thanasilp, J. Tangpanitanon, and D. G. Angelakis, Qubit-efficient encoding schemes for binary optimisation problems, Quantum 5, 454 (2021).
  * [49] E. X. Huber, B. Y. L. Tan, P. R. Griffin, and D. G. Angelakis, Exponential qubit reduction in optimization for financial transaction settlement, arXiv , 2307.07193.
  * [50] I. D. Leonidas, A. Dukakis, B. Y. L. Tan, and D. G. Angelakis, Qubit efficient quantum algorithms for the vehicle routing problem on nisq processors, arXiv , 2306.08507.
  * Azad  _et al._ [2022] U. Azad, B. K. Behera, E. A. Ahmed, P. K. Panigrahi, and A. Farouk, Solving vehicle routing problem using quantum approximate optimization algorithm, IEEE Trans. Intell. Transp. Syst. , 7564 (2022).
  * Perelshtein  _et al._ [2023] M. R. Perelshtein, A. I. Pakhomchik, A. A. Melnikov, M. Podobrii, A. Termanova, I. Kreidich, B. Nuriev, S. Iudin, C. W. Mansell, and V. M. Vinokur, Nisq-compatible approximate quantum algorithm for unconstrained and constrained discrete optimization, Quantum 7, 1186 (2023).
  * Lucas [2014] A. Lucas, Ising formulations of many np problems, Front. Phys. 2, 5 (2014).
  * [54] See supplementary materials.
  * Wales and Doye [1997] D. J. Wales and J. P. K. Doye, Global optimization by basin-hopping and the lowest energy structures of lennard-jones clusters containing up to 110 atoms, J. Phys. Chem. A 101, 5111 (1997).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, On the limited memory bfgs method for large scale optimization, Math. Program. 45, 503 (1989).
  * Glover  _et al._ [2010] F. Glover, Z. Lü, and J.-K. Hao, Diversification-driven tabu search for unconstrained binary quadratic problems, 4OR 8, 239 (2010).
  * Dunning  _et al._ [2018] I. Dunning, S. Gupta, and J. Silberholz, What works best when? a systematic evaluation of heuristics for max-cut and QUBO, INFORMS Journal on Computing 30, 608 (2018).
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices, Phys. Rev. X 10, 021067 (2020).
  * [60] J. Basso, E. Farhi, K. Marwaha, B. Villalonga, and L. Zhou, The quantum approximate optimization algorithm at high depth for maxcut on large-girth regular graphs and the sherrington-kirkpatrick model, arXiv , 2110.14206.
  * [61] A. Galda, E. Gupta, J. Falla, X. Liu, D. Lykov, Y. Alexeev, and I. Safro, Similarity-based parameter transferability in the quantum approximate optimization algorithm, arXiv , 2307.05420.
  * Brady  _et al._ [2021] L. T. Brady, C. L. Baldwin, A. Bapat, Y. Kharkov, and A. V. Gorshkov, Optimal protocols in quantum annealing and quantum approximate optimization algorithm problems, Phys. Rev. Lett. 126, 070505 (2021).
  * Chen  _et al._ [2022] Y. Chen, L. Zhu, N. J. Mayhall, E. Barnes, and S. E. Economou, How much entanglement do quantum optimization algorithms require?, in  _Quantum 2.0_ (Optica Publishing Group, 2022) pp. QM4A–2.
  * Dupont  _et al._ [2022] M. Dupont, N. Didier, M. J. Hodson, J. E. Moore, and M. J. Reagor, Entanglement perspective on the quantum approximate optimization algorithm, Phys. Rev. A 106, 022423 (2022).
  * Nakhl  _et al._ [2024] A. C. Nakhl, T. Quella, and M. Usman, Calibrating the role of entanglement in variational quantum circuits, Phys. Rev. A 109, 032413 (2024).
  * Parisi [1979] G. Parisi, Infinite number of order parameters for spin-glasses, Phys. Rev. Lett. 43, 1754 (1979).
  * Boettcher [2005] S. Boettcher, Extremal optimization for sherrington-kirkpatrick spin glasses, Eur. Phys. J. B 46, 501 (2005).

Supplementary Material for Qubit-efficient quantum combinatorial optimization
solver

##  S.1 Qubit-efficient encoding

In the main text, we described the qubit-efficient encoding of $N$ classical
variables [Eq. (1)], and how to estimate $\bar{z}_{i}$ and
$\overline{z_{i}z_{j}}$ [Eqs. (2) and (3)]. Here, we exemplify these using the
$(N=4,\ d=2)$ case as an example.

Encoding $N=4$ variables using $d=2$ requires 3 qubits. The most generic
3-qubit wave function is

$\ket{\psi}=c_{0}\ket{000}+c_{1}\ket{001}+\cdots+c_{7}\ket{111}$ (S.1)

with $\sum_{i=0}^{7}|c_{i}|^{2}=1$. Denoting the most significant bit as the
label bit, projectively measuring this wave function gives the following
measurements with respective probabilities,

$\begin{array}[]{c|c|c}\textrm{Variables}&\textrm{Values}&\textrm{Probability}\\\
\hline\cr(z_{0},z_{1})&(1,1)&|c_{0}|^{2}\\\
(z_{0},z_{1})&(1,-1)&|c_{1}|^{2}\\\ (z_{0},z_{1})&(-1,1)&|c_{2}|^{2}\\\
(z_{0},z_{1})&(-1,-1)&|c_{3}|^{2}\\\ (z_{2},z_{3})&(1,1)&|c_{4}|^{2}\\\
(z_{2},z_{3})&(1,-1)&|c_{5}|^{2}\\\ (z_{2},z_{3})&(-1,1)&|c_{6}|^{2}\\\
(z_{2},z_{3})&(-1,-1)&|c_{7}|^{2}\\\ \hline\cr\end{array}$ (S.2)

In an experiment with $N_{s}$ shots, $N_{s}^{(0)}\sim
N_{s}\sum_{i=0}^{3}|c_{i}|^{2}$ shots yield the values of $(z_{0},z_{1})$, and
$N_{s}^{(1)}\sim N_{s}\sum_{i=4}^{7}|c_{i}|^{2}$ shots yield the values of
$(z_{2},z_{3})$. The estimate for, e.g., $\overline{z_{0}z_{1}}$ if obtained
from just the $N_{s}^{(0)}$ shots as

$\displaystyle\overline{z_{0}z_{1}}=$ $\displaystyle\frac{\sum_{i,j=\pm
1}ijN(z_{0}=i,z_{1}=j)}{N_{s}^{(0)}}$ $\displaystyle\sim$
$\displaystyle\frac{N_{s}|c_{0}|^{2}-|c_{1}|^{2}-|c_{2}|^{2}+|c_{3}|^{2})}{N_{s}(|c_{0}|^{2}+|c_{1}|^{2}+|c_{2}|^{2}+|c_{3}|^{2})}$
(S.3)

where $N(z_{0}=i,z_{1}=j)$ is the number of times $(z_{0},z_{1})$ were
measured as $(i,j)$. The expression in Eq. (S.1) is equal to
$\frac{\braket{\hat{P}_{\ell=0}\hat{Z}_{1}\hat{Z}_{2}}}{\braket{\hat{P}_{\ell=0}}}$.
Other correlations and one-body terms can be similarly estimated as

$\displaystyle\overline{z_{0}z_{1}}=$
$\displaystyle\frac{|c_{0}|^{2}-|c_{1}|^{2}-|c_{2}|^{2}+|c_{3}|^{2}}{|c_{0}|^{2}+|c_{1}|^{2}+|c_{2}|^{2}+|c_{3}|^{2}}$
$\displaystyle\overline{z_{2}z_{3}}=$
$\displaystyle\frac{|c_{4}|^{2}-|c_{5}|^{2}-|c_{6}|^{2}+|c_{7}|^{2}}{|c_{4}|^{2}+|c_{5}|^{2}+|c_{6}|^{2}+|c_{7}|^{2}}$
$\displaystyle\bar{z}_{0}=$
$\displaystyle\frac{|c_{0}|^{2}+|c_{1}|^{2}-|c_{2}|^{2}-|c_{3}|^{2}}{|c_{0}|^{2}+|c_{1}|^{2}+|c_{2}|^{2}+|c_{3}|^{2}}$
$\displaystyle\bar{z}_{1}=$
$\displaystyle\frac{|c_{0}|^{2}-|c_{1}|^{2}+|c_{2}|^{2}-|c_{3}|^{2}}{|c_{0}|^{2}+|c_{1}|^{2}+|c_{2}|^{2}+|c_{3}|^{2}}$
$\displaystyle\bar{z}_{2}=$
$\displaystyle\frac{|c_{4}|^{2}+|c_{5}|^{2}-|c_{6}|^{2}-|c_{7}|^{2}}{|c_{4}|^{2}+|c_{5}|^{2}+|c_{6}|^{2}+|c_{7}|^{2}}$
$\displaystyle\bar{z}_{3}=$
$\displaystyle\frac{|c_{4}|^{2}-|c_{5}|^{2}+|c_{6}|^{2}-|c_{7}|^{2}}{|c_{4}|^{2}+|c_{5}|^{2}+|c_{6}|^{2}+|c_{7}|^{2}}$
(S.4)

The estimate for the cost $C=\sum_{ij}w_{ij}z_{i}z_{j}$ is

$\displaystyle C=$ $\displaystyle
w_{12}\overline{z_{1}z_{2}}+w_{34}\overline{z_{3}z_{4}}$
$\displaystyle+w_{13}\bar{z}_{1}\bar{z}_{3}+w_{14}\bar{z}_{1}\bar{z}_{4}+w_{23}\bar{z}_{2}\bar{z}_{3}+w_{24}\bar{z}_{2}\bar{z}_{4}.$
(S.5)

###  S.1.1 d=N

In the $d=N$ case, there are no label bits, and each variable is mapped one-
to-one with a data bit. This is the usual way variables are mapped to qubits.
Correlations $z_{i}z_{j}$ are, e.g., estimated as
$\overline{z_{i}z_{j}}=\braket{\hat{Z}_{i}\hat{Z}_{j}}$, and
$C=\sum_{ij}w_{ij}z_{i}z_{j}$ is $C=\braket{\hat{\mathcal{H}}}$ with
$\hat{\mathcal{H}}=\sum_{ij}w_{ij}\hat{Z}_{i}\hat{Z}_{j}$.

###  S.1.2 d=1

In the $d=1$ case, there are $\log_{2}N$ label bits, and all the variables are
mapped to one data bit. All correlations $z_{i}z_{j}$ are estimated as
$\overline{z_{i}z_{j}}=\bar{z}_{i}\times\bar{z}_{j}$.

###  S.1.3 d=2

The $d=2$ case, coincidentally, requires the same number of qubits as $d=1$.
However, the mappings between qubits and variables are different. For $d=2$,
there are $\log_{2}N-1$ label bits and $2$ data bits, unlike the $d=1$ case.

##  S.2 Entanglement in the qubit-efficient encoding

Entanglement is an essential property of quantum devices and algorithms, and
is the reason for non-simulability and potential utility of quantum
optimization. The growth of entanglement induced by the QAOA circuit has been
previously studied [63, 64, 65], finding for example volume-law entanglement
in the wave function between the initial and final state. The target state of
QAOA is a product state with no entanglement.

In the qubit-efficient approach, however, the target state is not a product
state. Instead, the bit string solution $\boldsymbol{z}$ is stored as the
entangled wave function
$\ket{\phi(\boldsymbol{z})}=\sum_{\ell}\lambda_{\ell}\ket{\ell}\otimes\ket{\boldsymbol{z}_{\ell}}$
[Eq. (1) in the main text]. $\boldsymbol{z}_{\ell}$ are length-$d$ subsets of
$\boldsymbol{z}$, where the subscript $\ell$ runs from $0$ to $N/d-1$. Here,
we quantify the entanglement present in the target wave function by
calculating the entanglement entropy of the data qubits.

The reduced density matrix for the data qubits, obtained by tracing over the
label qubits, is

$\rho_{\rm
data}=\sum\nolimits_{\ell=0}^{N/d-1}|\lambda_{\ell}|^{2}\ket{\boldsymbol{z}_{\ell}}\bra{\boldsymbol{z}_{\ell}},$
(S.6)

and its von Neumann entropy is $S_{\rm data}=-{\rm Tr}(\rho_{\rm
data}\log_{2}\rho_{\rm data})$.

Figure S.1: Average entanglement entropy of the data qubits that store bit
strings of length $N=2^{16}$. The data qubits have volume law entanglement.

Figure S.1 plots the average entanglement entropy, averaged over 100 random
bit strings of length $N=2^{16}$, versus $d$, assuming wave functions with
uniform $|\lambda_{\ell}|^{2}=d/N$. Two regimes emerge. When $d\lesssim 16$,
the entanglement entropy increases as $d$, whereas for $d\gtrsim 16$, it
decreases as $\log_{2}(N/d)$. This entropy is the smaller of the number of
data and label qubits, indicating volume law entanglement.

We can understand the two regimes as follows. When we encode a random bit
string $\boldsymbol{z}$ at large $d$, i.e. $2^{d}\gg N/d$, the probability
that any two $\boldsymbol{z}_{\ell}$ are equal is exponentially small in $d$,
therefore on average, $S_{\rm
data}=-\sum_{\ell=0}^{N/d-1}|\lambda_{\ell}|^{2}\log_{2}|\lambda_{\ell}|^{2}$.
For uniform $|\lambda_{\ell}|^{2}=d/N$, the entropy reduces to $S_{\rm
data}=\log_{2}(N/d)$.

When $2^{d}\ll N/d$, there is a significant probability that there will be
repetitions among the different $\boldsymbol{z}_{\ell}$. The average number of
times any $\boldsymbol{z}_{\ell}$ would be repeated in $\boldsymbol{z}$ is
$\frac{N/d}{2^{d}}$. Those repetitions must be grouped together in $\rho_{\rm
data}$. For uniform $\lambda_{\ell}^{2}=d/N$, this gives $\rho_{\rm
data}\sim\sum_{\boldsymbol{z}_{\ell}=0}^{2^{d}-1}\frac{d}{N}\frac{N/d}{2^{d}}\ket{\boldsymbol{z}_{\ell}}\bra{\boldsymbol{z}_{\ell}}$.
The entropy of this state is $S_{\rm data}=d$.

##  S.3 Example problem in Fig. 2

In the main text, we exemplified the qubit-efficient encoding [Eq. (1)] and
the ansatz using an SK problem with $N=4$ variables. Here, we discuss this
example in detail.

###  S.3.1 Example problem

The specific problem considered in Fig. 1 in the main text has the weights

$w=\left(\begin{array}[]{cccc}0&1&-1&1\\\ 1&0&-1&-1\\\ -1&-1&0&1\\\
1&-1&1&0\end{array}\right).$ (S.7)

It can be verified, e.g. by explicitly enumerating all the 16 bitstrings for
the 4 variables, that the ground state solutions for this example are
$\boldsymbol{z}^{*}=(1,-1,1,-1)$ and $-\boldsymbol{z}^{*}=(-1,1,-1,1)$, with
the ground state cost $C^{*}=-4$. These solutions can be represented by the
respective wave functions
$\ket{\phi(\boldsymbol{z}^{*})}=(\ket{0}_{\textrm{label}}\otimes\ket{01}_{\textrm{data}}+\ket{1}_{\textrm{label}}\otimes\ket{01}_{\textrm{data}})/\sqrt{2}$
and
$\ket{\phi(-\boldsymbol{z}^{*})}=(\ket{0}_{\textrm{label}}\otimes\ket{10}_{\textrm{data}}+\ket{1}_{\textrm{label}}\otimes\ket{10}_{\textrm{data}})/\sqrt{2}$
at $d=2$.

To implement the ansatz
$\ket{\psi_{p=1}}=e^{-i\beta\hat{\mathcal{H}}_{x}/2}e^{-i\gamma\hat{\mathcal{H}}[\psi_{0}]/2}\ket{\psi_{0}}$,
we first need to determine $\hat{\mathcal{H}}[\psi_{0}]$, where
$\ket{\psi_{0}}$ is the equal superposition state of 3 qubits. Measurements of
$\ket{\psi_{0}}$ would give

$\displaystyle\bar{z}_{1}=\bar{z}_{2}=\bar{z}_{3}=\bar{z}_{4}=0,$
$\displaystyle\braket{\psi_{0}}{\hat{P}_{\ell=0}\psi_{0}}=\braket{\psi_{0}}{\hat{P}_{\ell=1}\psi_{0}}=1/2,$
$\displaystyle\Rightarrow$
$\displaystyle\hat{\mathcal{H}}[\psi_{0}]=2(w_{12}\hat{P}_{\ell=0}\hat{Z}_{d_{1}}\hat{Z}_{d_{2}}+w_{34}\hat{P}_{\ell=1}\hat{Z}_{d_{1}}\hat{Z}_{d_{2}}).$
(S.8)

Next, we describe how to implement the ansatz at arbitrary $p$, which can be
straightforwardly specialized to $p=1$.

###  S.3.2 Circuit compilation of the ansatz

Figure S.2: Simplification of the phase separator term in the ansatz.
Implementation of these logical gates in terms of hardware-native gates are
given in Fig. S.3

Here, we describe how to implement the ansatz
$\ket{\psi_{p}}=e^{i\beta\hat{\mathcal{H}}_{x}}e^{i\gamma\hat{\mathcal{H}}[\psi_{p-1}]}\ket{\psi_{p-1}}$
for the specific case of $N=4,d=2$ considered here. A general implementation
of the ansatz for arbitrary $N$ and $d$ is discussed later.

Figure S.2 breaks down $e^{i\gamma\hat{\mathcal{H}}[\psi_{p-1}]/2}$ into a
3-qubit gate and 2-qubit gate. Figure S.3 shows how to implement compile these
gates into the native gates on Rigetti’s AnkaaTM-9Q-3 chip,
$R_{x}(\pi/2)\equiv\exp(-i\hat{X}\pi/4)$ and
$\textrm{iSWAP}=\exp(i\pi/4(\hat{X}\otimes\hat{X}+\hat{Y}\otimes\hat{Y}))$.
$R_{z}(\phi)$ gates are virtually applied via frame tracking. The unitary
$e^{i\beta\hat{\mathcal{H}}_{x}}$ can be implemented with $R_{x}(-2\beta)$ on
all the qubits. In Fig. S.3,

$\displaystyle U_{1}=R_{x}(-\pi/2)$ $\displaystyle
U_{2}(\phi)=R_{z}(\pi/2)R_{x}(\pi/2)R_{z}(\phi)R_{x}(-\pi/2)$ $\displaystyle
U_{3}=R_{z}((\pi+\phi)/2)R_{x}(\pi/2)$ $\displaystyle
U_{4}=R_{z}(\pi/2)R_{x}(\pi/2)$ (S.9)

Figure S.3: Compilation of common gates in terms of gates native to Rigetti’s
AnkaaTM-9Q-3.

###  S.3.3 Optimization of the ansatz at p=1

Figure S.4: Distribution of bit strings measured on Rigetti’s AnkaaTM-9Q-3
chip (teal), and predicted from the wave function (magenta), at the optimal
angles of the ansatz for the example in Fig. 2 in the main text.

After implementing the ansatz
$\ket{\psi_{p=1}}=e^{i\beta_{1}\hat{\mathcal{H}}_{x}}e^{i\gamma_{1}\hat{\mathcal{H}}[\psi_{0}]}\ket{\psi_{0}}$,
explicit calculation of the expectation value of the cost gives

$\displaystyle C[\psi_{p=1}]=$ $\displaystyle(w_{12}-w_{34})\sin
2\beta_{1}(\sin(4\gamma_{1}w_{12})-\sin(4w_{34}\gamma_{1})+\sin(2\gamma_{1}(w_{12}-w_{34})))$
$\displaystyle+\sin 2\beta_{1}\cos
2\beta_{1}(w_{12}+w_{34})(\sin(4w_{12}\gamma_{1})+\sin(4w_{34}\gamma_{1})).$
(S.10)

For the specific example in Eq. (S.7), the noiseless cost is
$C[\psi_{p=1}]=2\sin 4\beta_{1}\sin 4\gamma_{1}$. The cost measured on
AnkaaTM-9Q-3 as well as numerically computed from 10000 shots of a noiseless
simulation, plotted in Figs. 2(b-c) in the main text, closely match this
expression, with the former having an overall reduction in magnitude
attributable to hardware noise.

The parameters which minimize the cost at $p=1$ are
$(\beta_{1},\gamma_{1})=\left(\frac{3\pi}{8},\frac{\pi}{8}\right)$, and the
wave function $\ket{\psi_{p=1}}$ is

$\displaystyle\ket{\psi_{p=1}}=$
$\displaystyle\frac{(\ket{0}_{\textrm{label}}+\ket{1}_{\textrm{label}})\otimes(\ket{01}_{\textrm{data}}+\ket{10}_{\textrm{data}})}{2}$
$\displaystyle=$
$\displaystyle\frac{\ket{\phi(\boldsymbol{z}^{*})}+\ket{\phi(-\boldsymbol{z}^{*})}}{\sqrt{2}}.$
(S.11)

The wave function is an equal superposition of
$\ket{\phi(\boldsymbol{z}^{*})}$ and $\ket{\phi(-\boldsymbol{z}^{*})}$ due to
$\mathbb{Z}_{2}$ symmetry of the ansatz. Figure S.4 plots the distribution of
bit strings sampled from this wave function (magenta), and measured in
executing the ansatz on Rigetti’s AnkaaTM-9Q-3 chip. The ideal wave function
gives

$\displaystyle\overline{z_{1}z_{2}}=\overline{z_{3}z_{4}}=-1,$
$\displaystyle\bar{z}_{1}=\bar{z}_{2}=\bar{z}_{3}=\bar{z}_{4}=0,$ (S.12)

therefore the cost is $C[\psi_{p=1}]=-(w_{12}+w_{34})=-2$. The measured cost
on Rigetti’s AnkaaTM-9Q-3 is $C=-1.736(1)$. We remark on the unique feature of
the qubit-efficient encoding that a superposition of ground state solutions
does not have the same cost $C$ as the ground state cost $C^{*}$, whereas the
cost for the superposition would be equal to the ground state cost $C^{*}$ in
the $d=N$ case. This is because $\bar{z}_{i}=0$, therefore several terms in
the cost at $d<N$ are 0. Therefore, explicitly breaking the $\mathbb{Z}_{2}$
symmetry in the ansatz is necessary. We broke the symmetry with the one-body
bias term, $\exp(i\gamma^{\prime}\hat{\mathcal{H}}_{z})$.

##  S.4 Parameter clustering in the SK model and the qubit-efficient ansatz

Figure S.5: Optimal parameters returned by the classical optimizer for the
five SK problem instances with $N=64$ variables in Fig. 3 in the main text, at
$d=4$. The thin curves plot the instance-by-instance optimal parameters, and
the thick curves plot their mean values. We find the curves are clustered
together.

Previous works [43, 46] have found evidence of parameter clustering when
solving SK problem instances using QAOA (which is the $d=N$ limit of our
ansatz). Optimal parameters and approximation ratios for the typical problem
instance have been tabulated up to $p=12$. Our cost function $C$ is related to
the cost $\tilde{C}$ in Ref. [43] via $C=\sqrt{N}\tilde{C}$, and our
parameters $(\beta,\gamma)$ at $d=N$ are related to parameters
$(\tilde{\beta},\tilde{\gamma})$ in Ref. [43] via $\beta=-\tilde{\beta}$ and
$\gamma=-\tilde{\gamma}/\sqrt{N}$. The minus sign is irrelevant since $C$ is
real. The optimal value of $\gamma$, in our convention, scales as
$1/\sqrt{N}$. The optimal values of $(\beta,\gamma)$ and $C$ up to $p=11$ are
the same for two classes of SK models – $w_{ij}\in\\{\pm 1\\}$, and
$w_{ij}\sim\mathcal{N}(0,1)$.

We showed numerical evidence in Fig. 4 in the main text that parameter
clustering occurs for the qubit-efficient encoding ansatz ($d<N$) as well.
Figure S.5 plots the optimal parameters found by the classical optimizer for
the 5 problem instances in Fig. 3, up to $p=3$ at $d=4$ and $N=64$. While Ref.
[43] did not have – nor require (see previous section) – $\gamma^{\prime}$, we
introduced one-body bias terms $\gamma^{\prime}\hat{\mathcal{H}}_{z}$ to break
$\mathbb{Z}_{2}$ symmetry. We find that the thin curves, which plot the
optimal parameters for the 5 instances, are clustered together. We also find
that the optimal values of $\gamma$ are smaller than those of $\beta$ and
$\gamma^{\prime}$.

We hypothesize that the optimal coefficient multiplying
$\hat{Z}_{i}\hat{Z}_{j}$ in the qubit-efficient ansatz scales $\propto
1/\sqrt{N}$, just like in the $d=N$ case. We validate this hypothesis as
follows. Ising terms such as $\hat{Z}_{i}\hat{Z}_{j}$ appear in the ansatz via
$\gamma\frac{\hat{P}_{\ell}\hat{Z}_{d_{i}}\hat{Z}_{d_{j}}}{\braket{\hat{P}_{\ell}}}$.
Here, $\braket{\hat{P}_{\ell}}=d/N$ in the equal superposition state which is
the $p=0$ state of the ansatz, although it will vary with $p$. Setting
$\gamma\frac{N}{d}\propto 1/\sqrt{N}$ gives $\gamma(N,d)\propto 1/N^{3/2}$.
This is the rescaling we applied when we used optimal parameters from $N=64$
to solve SK problems instances with $N=128$ in Fig. 3(c), i.e.
$\gamma^{*}(N=128,d=4)=\gamma^{*}(N=64,d=4)/2^{3/2}$. We perform the same
rescaling and calculate the approximation ratios for an ensemble of 100 SK
problem instances with $N=256,512$ in Fig. S.6, and find the approximation
ratio increases with $p$.

Further, we validate the scaling of $\gamma$ with the following optimization
procedure. For an ensemble of SK problems at different sizes and $d$, we fix
the values of $\beta(N,d;p,i)$ and $\gamma^{\prime}(N,d;p,i)$ from one
instance at $N=64,d=4$, and vary the values of $\gamma(N,d;p,i)$, where $1\leq
i\leq p$. We assume that $\gamma(N,d;p,i)/\gamma(N=64,d=4;p,i)=\theta(N,d;p)$
is a constant, and optimize over $\theta(N,d;p)$. Figure S.7 finds that
$\theta^{\rm opt}(N,d;p)\propto d/N^{3/2}$, as predicted by the above
argument.

Figure S.6: Average approximation ratios for optimized parameters for 100 SK
problem instances at different sizes and $d=4$, where the same parameters were
used for all the instances. Figure S.7: Optimal ratios $\gamma^{\rm
opt}(N,d;p,i)/\gamma(N=64,d=64;p,i)$ scale as $d/N^{3/2}$. The points include
$1\leq p\leq 3$, $N\in\\{64,128,256\\},d\in\\{1,2,4,8\\}$, with one SK problem
instance at each $N$ and $d$. Dashed line plots $d/N^{3/2}$ as guide to the
eye.

##  S.5 Baseline performance of a related ansatz

Here, we consider a slightly modified ansatz given in Fig. S.8(a) The
difference with the ansatz in the main text is that there are no $R_{x}$
rotations of the label bits. Figure S.8(a) can be simplified to the ansatz in
Fig. S.8(b). This simplification uses the identity that applying a controlled
gate – apply $U$ on target bits if control bits’ value is $x$ – followed by a
measurement of the control bits is equivalent to first doing a measurement on
the control bits and then implementing $U$ on the target bits if the
measurement is $x$. The ansatz in Fig. S.8(b) essentially implements the QAOA
ansatz on $d$ data bits using the SK problem weights on the subset of the $d$
variables corresponding to the measured label value. Repeating the ansatz in
Fig. S.8(b) several times essentially solves all the $d$-variable subproblems
with label $\ell$, obtains the cost $C_{\ell}$ for the subproblems
$\ell=0\cdots(N/d-1)$, and returns the final cost,

$C=\sum_{\ell=0}^{N/d-1}C_{\ell}.$ (S.13)

Figure S.8: An ansatz with phase separators and one-qubit rotations on data
qubits only. Executing the ansatz in (a) is equivalent to repeating the ansatz
in (b) for $\ell=0\cdots(N/d-1)$, and reduces to implementing QAOA on the
groups with $d$ variables each.

When the size of the problem being solved by QAOA is large, Ref. [43] predicts
clustering of QAOA parameters and tabulates the approximation ratio,
$r^{*}(p)$ up to $p=11$. The asymptotic value of the ground state cost for $N$
variables is $C^{*}_{\textrm{asymp.}}(N)=-PN^{3/2}+\alpha N^{\omega}$, where
$P\approx 0.7632$ is the Parisi constant [66], $\omega=5/6$ and $\alpha\approx
0.7$ is a non-universal constant that accounts for finite-size effects [67].
Similarly, for $d\gg 1$, the ground state cost for the subproblems with $d$
variables is $C^{*}_{\textrm{asymp.}}(d)\sim-Pd^{3/2}$, neglecting the
subleading term, and the ansatz in Fig. S.8 at $p$ layers yields
$C_{\ell}^{*}(p)=r^{*}(p)\times Pd^{3/2}$. Then, the total cost for $N$
variables is

$C^{*}=r^{*}(p)\times Pd^{3/2}\times\frac{N}{d}.$ (S.14)

The approximation ratio for this ansatz is

$r^{*}\equiv C^{*}/C^{*}_{\textrm{asymp.}}(N)\sim
r^{*}(p)\times\sqrt{\frac{d}{N}}.$ (S.15)

This approximation ratio is plotted as the black dash-dotted lines in Fig. 4.

##  S.6 Circuit implementation of the ansatz

Figure S.9: The phase separator unitary in the ansatz [Eq.(6) in the main
text] is equivalent to implementing controlled unitaries with control as label
qubits, and targets as Ising interactions on the data qubits. Figure S.10:
Reduction of controlled-unitary with $m$ control bits to a controlled-unitary
with $(m-1)$ control bits and one ancilla bit. This simplification, along with
the compilation of the Toffoli gate in (b), can be recursively used to reduce
the controlled-unitary in Fig. S.9 to hardware-native gates.

###  S.6.1 Phase separator unitary

The phase separator Hamiltonian [Eq.(5)] in the main text is

$\displaystyle\hat{\mathcal{H}}[\psi]=$
$\displaystyle\sum_{\begin{subarray}{c}i<j\\\
\ell_{i}=\ell_{j}\end{subarray}}w_{ij}\frac{\hat{P}_{\ell_{i}}\hat{Z}_{d_{i}}\hat{Z}_{d_{j}}}{\braket{\psi}{\hat{P}_{\ell_{i}}}{\psi}}$
$\displaystyle+$ $\displaystyle\frac{1}{2}\sum_{\begin{subarray}{c}i<j\\\
\ell_{i}=\ell_{j}\end{subarray}}w_{ij}\left(\bar{z}_{i}\frac{\hat{P}_{\ell_{j}}\hat{Z}_{d_{j}}}{\braket{\psi}{\hat{P}_{\ell_{j}}}{\psi}}+\bar{z}_{j}\frac{\hat{P}_{\ell_{i}}\hat{Z}_{d_{i}}}{\braket{\psi}{\hat{P}_{\ell_{i}}}{\psi}}\right).$
(S.16)

Since all the terms in $\hat{\mathcal{H}}[\psi]$ commute with each other, the
unitary $\exp(i\gamma\hat{\mathcal{H}}[\psi])$ is a product of a unitary due
to each term. The unitary due to each term is of the form
$\exp(i\gamma\hat{P}_{\ell}\hat{O}_{\textrm{d}ata})$ where
$\hat{O}_{\textrm{data}}=\hat{Z}_{d_{i}}$ or
$\hat{O}_{\textrm{data}}=\hat{Z}_{d_{i}}\hat{Z}_{d_{j}}$. When the label
qubits are in the state $\ket{\ell}$, the projection operator $\hat{P}_{\ell}$
is 1, therefore $\exp(i\gamma\hat{P}_{\ell}\hat{O}_{\textrm{d}ata})$
implements $\exp(i\gamma\hat{O}_{\textrm{d}ata})$ on the data qubits. When the
label qubits are in $\ket{\ell^{\prime}\neq\ell}$, the projection operator
$\hat{P}_{\ell}$ is 0, therefore
$\exp(i\gamma\hat{P}_{\ell}\hat{O}_{\textrm{d}ata})$ does nothing. This is
nothing but a controlled-unitary operator, which implements the target unitary
$\exp(i\gamma\hat{O}_{\textrm{d}ata})$ if the label qubits are in
$\ket{\ell}$, as depicted in Fig. S.9. One method to implement controlled
unitaries with many control bits is shown in Fig. S.10(a).

##  S.7 Shot noise

In Fig. 3 of the main text, we variationally optimized the cost computed from
a noiseless calculation of the wave function using a state vector simulator.
In Figs. S.11 and S.12, we show that the estimated cost converges to the exact
cost from the wave function at $O(1000)$ shots for $N=64$ and $128$.

Figure S.11: Estimated cost versus number of shots for the 5 SK problem
instances in Fig. 3 of the main text at $N=64,d=4$. The costs converge to the
exact cost computed from the wave function (dashed line) at $\sim 1000$ shots.
Figure S.12: Estimated cost versus number of shots for the 5 SK problem
instances in Fig. 3 of the main text at $N=128,d=4$. The costs converge to the
exact cost computed from the wave function (dashed line) at $O(1000)$ shots.