Fast Beam Training and Performance Analysis for Extremely Large Aperture Array This work was supported by the National Key R&D Project of China under Grant No. 2023YFB2905000. (Corresponding Author: Hongkang Yu)

Yuan Si12,
Hongkang Yu12,
Yijian Chen12,

1State Key Laboratory of Mobile Network and Mobile Multimedia Technology, Shenzhen 518055, China

2Wireless and Computing Product R&D Institute, ZTE Corporation, Xian 710114,  China

{si.yuan, yu.hongkang<EMAIL_ADDRESS>
Extremely large aperture array (ELAA) can significantly enhance beamforming gain and spectral efficiency. Unfortunately, the use of narrower beams for data transmission results in a substantial increase in the cost of beam training. In this paper, we study a high-efficiency and low-overhead scheme named hash beam training. Specifically, two improved hash codebook design methods, random and fixed, are proposed. Moreover, we analyze beam alignment performance. Since the derived beam alignment success probability is a complex function, we also propose a heuristic metric to evaluate the impact of codebook parameter on performance. Finally, simulation results validate the theoretical analysis, indicating that the proposed beam training scheme can achieve fast beam alignment with lower overhead and higher accuracy.

ELAA, beam training, codebook design, hash function, performance analysis.

§ INTRODUCTION

With the development of wireless communication, the operating frequency becomes higher, and the array aperture becomes larger. 5G new radio (NR), which uses massive multi-input multi-output (MIMO) as one of the key technologies, has applied the millimeter-wave band and supports 32 antenna ports with hundreds of antennas [1]. In the future, 6G will further introduce terahertz (THz) wave band and the extremely large aperture array (ELAA) equipped with thousands of antennas. More antennas and larger array aperture can bring higher beamforming gain to combat the severe path loss in high frequency bands and ensure the coverage of the base station [2].

However, higher beam gain leads to narrower beams, which makes the overhead of beam training overwhelming. How to achieve fast beam alignment with lower overhead is a major challenge. The traditional beam sweeping scheme sends training beams sequentially, which brings high overhead, is not suitable for the ELAA scenario. To solve this, hierarchical beam training schemes were proposed [3, 4, 5], where the base station first sends few wide beams, narrows the beam search range according to the feedback from the user, and finally achieves beam alignment with specific resolution. Unfortunately, this scheme can reduce overhead only in the case of few users, and requires multiple feedbacks, which increases the delay.

Channel estimation based on compressive sensing is another scheme that can achieve high-precision beam alignment [6, 7]. Taking advantage of the sparsity  of the channel in the angular domain at high frequency bands, we can obtain the estimated channel vector by the sparse signal recovery algorithm with a small amount of pilot overhead. However, the low overhead comes at the expense of robustness and increased computational complexity. Moreover, the channel estimation performance deteriorates rapidly as the signal-to-noise ratio (SNR) decreases.

Recently, a high-efficiency scheme named hash beam training received extensive attention, where the training beams comprise multiple main lobes with different directions [8, 9]. Based on the measurement results, the user can infer whether the optimal beam direction is among these specific main lobes. Through cross-validation of multiple measurement results, the best beam direction can be selected and fed back to the base station. Hash beam training not only saves training overhead but also has lower computational complexity, enabling better support for multi-user scenarios and applicable in both ELAA and reconfigurable intelligent surfaces (RIS). The design of the training beam codebook has a significant impact on the performance of the hash beam training scheme. Unfortunately, existing hash codebooks are generated randomly without specialized design, so there is still a large space for improvement.

Addressing the aforementioned issues, this paper focuses on the codebook design in hash beam training, proposing two improved codebook design schemes, random and fixed, and conducts a theoretical analysis of the beam alignment success probability. Due to the complex form of the derived theoretical results, a novel heuristic metric with closed-form expressions is also proposed to evaluate the impact of hash codebooks on beam alignment performance. According to this parameter, the key parameters used to generate the hash codebook can be quickly determined. Finally, simulation results verify the correctness of the theoretical analysis and show better performance than the existing schemes.

§ SYSTEM MODEL

Consider a THz point-to-point multi-input single-output (MISO) system, where the base station is equipped with a uniform linear array (ULA) of $N$ antennas to serve a single-antenna user. Assuming that the user position is constant or keep moving at a low speed, the channel has the block fading property, that is, it can be regarded as constant over a period of time. Considering the sparsity of the THz channel, we use the geometric channel model [2], which can be expressed by

\begin{equation}
\bold{h}=\sum^P_{p=1} \alpha_p \bold{a}(\theta_p),
\end{equation}

where $\bold{a}(\theta)=\frac{1}{\sqrt{N}}[1,e^{-j\pi\sin\theta},\dots,e^{-j\pi(N-1)\sin\theta}]^{T}\in\mathbb{C}^{N\times1}$  represents the array response vector, $P$ denotes the number of paths, $\alpha_p$ and $\theta_p$ are the gain and angle of the $p$-th path, respectively. For THz band, the value of $P$ is usually small, e.g., 2 to 3 [2].

We assume that the purpose of beam training is to select an optimal beam from $N$ discrete Fourier transform (DFT) beams for data transmission, the beamforming weight of the $n$-th $(n=1,2,\dots,N)$ DFT beam can be expressed as

\begin{equation}
\bold{f}_n=\bold{a}(\theta_n),
\end{equation}

where $\mathrm{sin}\theta_n=-1+2n/N$.

In hash beam training scheme, the base station sequentially sends $M$ training beams, each of which is a superposition of $L$ DFT beams, as shown in Fig. <ref>. The indices of the selected $L$ DFT beams are randomly generated by a hash function, and we use the hash codebook $\mathbf{C}\in\mathbb{C}^{M\times N}$ to describe this process. If $c_{m,n}=1$, it means that the $n$-th DFT beam $\bold{f}_n$ is selected by the $m$-th training beam, while $c_{m,n}=0$ denotes that the $n$-th DFT beam is not selected. Each training beam selects $L$ DFT beams, i.e., $\sum^N_{n=1}c_{m,n}=L$, and the beamforming weight of the $m$-th training beam can be expressed as

\begin{equation}
\bold{w}_m=\frac{1}{\sqrt{L}}\sum^N_{n=1}c_{m,n}\bold{f}_n.
\end{equation}

Based on the above model, the reference signal received power (RSRP) of the user in the $m$-th training can be expressed as

\begin{equation}\label{eq4}
y_m=\left| \bold{h}^H\bold{w}_m \right|^2+n_m,
\end{equation}

where $n_m\sim \mathcal{N}(0,\sigma^2)$ denotes Gaussian noise.

The advantage of hash-based beam training is that the user can identify the best DFT beam based only on the RSRP $\{y_1,y_2,\dots,y_M\}$, which leads to the low processing complexity. Since the same hash function is preset in the user and the base station, the user can obtain the hash codebook $\bold{C}$ according to the training beam index $m$. After that, users select the best beam by voting mechanism. If the $m$-th training contains the $n$-th DFT beam, the received power $y_m$ of the $m$-th training is accumulated to the measured result $p_n$ of the DFT beam $n$, then the measurement result of the $n$-th beam $p_n$ can be expressed as

\begin{equation}\label{eq5}
\end{equation}

Finally, the user selects the DFT beam $\bold{f}_{\hat{n}}$ with the largest measurement result for reporting, where $\hat{n}=\arg\mathop{\max}\limits_{n} p_n$.

Comment 1: Many existing beam training schemes can be regarded as special cases of hash beam training, and the difference is only in the hash codebook. For example, the beam sweeping scheme can be represented by a hash codebook $\bold{C}=\bold{I}$ with $M=N,L=1$. The hierarchical beam training scheme in [10] corresponds to a hash codebook with $M=2\log_{2}N,L=N/2$. Therefore, the hash codebook design has a significant impact on the beam alignment performance as well as the training overhead, which will be analyzed in the following two sections.

An illustration of hash beam training scheme.

§ HASH CODEBOOK DESIGN

This section proposes several improved codebook design methods, incorporating both random and fixed approaches. In the existing scheme described in [8], each training beam independently selects DFT beams. Given the limited number of training beams, there is a significant variance in the selection frequency of individual DFT beams ${G_n} = \sum\nolimits_{m = 1}^M {{c_{m,n}}} $. As a result, when users are located in different positions, the probability of successful selection for each DFT beam varies considerably, leading to a decrease in overall beam alignment performance.

§.§ The Proposed Random Hash Codebook

Addressing the issues present in the existing scheme, the proposed random hash codebook selects DFT beams through a collaborative approach for each training beam, with the goal of equalizing the selection frequency $G_n$ of each DFT beam. To achieve this objective, we constrain the codebook parameters such that $N$ is divisible by $L$, and $M$ is divisible by $N/L$. We divide the training beams into $ML/N$ groups, each containing $N/L$ training beams. Within the same group, training beams randomly select DFT beams through a sampling without replacement, ensuring that each DFT beam appears only once in a training beam group, thereby satisfying ${G_n} = ML/N,\forall n$. The detailed process of the proposed scheme is summarized in Algorithm <ref>.
Proposed Random Hash Codebook

[1]

Number of DFT beams $N$, number of training beams $M$, codebook parameter $L$ Hash codebook $\mathbf C$
Initialize ${\mathbf{C}} = {{\mathbf{0}}_{M \times N}}$
$\bmod \left( {m - 1,N/L} \right) =  = 0$
$\mathcal{S} = \left\{ {1,2, \dots ,N} \right\}$
Randomly select $L$ elements from $\mathcal S$ to form set $\mathcal{T}$
Set ${c_{m,n}} = 1$ for $n\in\mathcal T$
$\mathcal{S} = \mathcal{S} - \mathcal{T}$

§.§ The Proposed Fixed Hash Codebook

When considering system implementation, the use of random hash codebooks leads to increased complexity. For example, base stations and users need to pre-configure the same random seed for codebook generation. To address this, this subsection introduces the design of a fixed hash codebook.

Since the codebook is predetermined in actual scenarios, it is feasible to search for the optimal codebook for beam alignment performance using offline simulation. Specifically, we determine the user distribution based on actual scenarios to obtain channel samples, and utilize the method described in Section <ref> to generate $X$ proposed random hash codebooks $\left\{ {{{\mathbf{C}}_1},{{\mathbf{C}}_2}, \ldots ,{{\mathbf{C}}_X}} \right\}$. Subsequently, we perform Monte Carlo simulations on these $X$ hash codebooks under identical channel samples, and select the one with the highest beam alignment accuracy as the fixed hash codebook.

Compared to randomly generated hash codebooks, the proposed fixed hash codebook is not only easier to implement in the system but also better matches the actual environment. Related simulation results will be presented in Section <ref>.

\begin{equation}\label{top1}
\begin{split}
\mathbb{P}\left\{ {{G_{{n^{*}}}} = {g_{{n^{*}}}},{G_{{n^{'}}}} = {g_{{n^{'}}}},{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}},} \right\} &= \\ \mathbb{P}\left\{ {{G_{{n^{*}}}} = {g_{{n^{*}}}}} \right\}\mathbb{P}\left\{ {{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}}|{G_{{n^{*}}}} = {g_{{n^{*}}}}} \right\}\mathbb{P}&\left\{ {{G_{{n^{'}}}} = {g_{{n^{'}}}}|{G_{{n^{*}}}} = {g_{{n^{*}}}},{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}}} \right\}
\end{split}
\end{equation}

\begin{equation}\label{final_result1}
\begin{split}
\tilde p = \sum\limits_{{g_{{n^{*}}}} = 0}^M {\sum\limits_{{k_{{n^{*}},{n^{'}}}} = 0}^{{g_{{n^{*}}}}} {\sum\limits_{{g_{{n^{'}}}} = {k_{{n^{*}},{n^{'}}}}}^M {\mathbb{P}\left\{ {{p_{{n^{*}}}} > {p_{{n^{'}}}}|{\mathbf{C}}} \right\}\mathbb{P}\left\{ {{G_{{n^{*}}}} = {g_{{n^{*}}}},{G_{{n^{'}}}} = {g_{{n^{'}}}},{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}},} \right\}} } }
\end{split}
\end{equation}

\begin{equation}\label{top2}
\begin{split}
\mathbb{P}\left\{ {{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}}} \right\} =  C_G^{{k_{{n^{*}},{n^{'}}}}}\left( {1 - \frac{{ML - NG}}{N}} \right)\left( {1 - \frac{{ML - NG}}{{N - 1}}} \right){\left( {\frac{{L - 1}}{{N - 1}}} \right)^{{k_{{n^{*}},{n^{'}}}}}}{\left( {1 - \frac{{L - 1}}{{N - 1}}} \right)^{G - {k_{{n^{*}},{n^{'}}}}}}, \\
{k_{{n^{*}},{n^{'}}}} = 0,1, \ldots ,G.
\end{split}
\end{equation}

§ BEAM ALIGNMENT PERFORMANCE ANALYSIS

This section analyzes the impact of hash codebooks on the beam alignment performance from a theoretical perspective. For the sake of simplicity, it is assumed that the channel contains only line-of-sight (LoS) components, i.e., $P=1$. This assumption is reasonable for high-frequency channels, as the gain of non-LoS (NLoS) paths is much smaller than that of LoS paths, having a minor impact on the received power. Furthermore, without loss of generality, we assume that the path gain satisfies $\alpha=1$.

§.§ Accurate Result of Beam Alignment Success Probability

Denoting the index of the optimal DFT beam as ${n^{*}}$, the probability of successful beam alignment can be represented as

\begin{equation}
\mathbb{P}\left\{ {{n^{*}} = \hat n} \right\} = \mathbb{P}\left\{ {\mathop  \cap \limits_{{n^{'}} \ne {n^{*}}}  {{p_{n^*}} > {p_{{n^{'}}}}}} \right\},
\end{equation}

where we focus on analyzing the probability of ${p_{{n^{*}}}} > {p_{{n^{'}}}}$. For the LoS channel ${\mathbf{h}} = {\mathbf{a}}\left( \theta  \right)$, ${\left| {{{\mathbf{h}}^{\text{H}}}{{\mathbf{f}}_n}} \right|^2}$ assumes a larger value only when $\theta$ is within the coverage range of the DFT beam. To simplify the derivation, it is further assumed that $\theta$ is located on a discrete grid point, i.e., $\theta  =  - 1 + 2n/N,n = 1,2, \dots ,N$, in which case ${n^{*}} = n$, and ${\left| {{{\mathbf{h}}^{\text{H}}}{{\mathbf{w}}_m}} \right|^2} = {c_{m,{n^{*}}}}/L$. Combining (<ref>) and (<ref>), we have

\begin{eqnarray}\label{eq7}
\begin{aligned}
{p_{{n^*}}} =& \frac{{{G_{{n^{*}}}}}}{L} + \sum\limits_{m = 1}^M {{c_{m,{n^{*}}}}{n_m}}, \\
{p_{{n^{'}}}} =& \frac{{{K_{{n^{*}},{n^{'}}}}}}{L} + \sum\limits_{m = 1}^M {{c_{m,{n^{'}}}}{n_m}} ,\forall {n^{'}} \ne {n^{*}},
\end{aligned}
\end{eqnarray}

where ${K_{{n^{*}},{n^{'}}}} = \sum\nolimits_{m = 1}^M {{c_{m,{n^{*}}}}{c_{m,{n^{'}}}}}$ represents the number of times DFT beams ${n^{*}}$ and ${n^{'}}$ are simultaneously selected by the training beams. According to (<ref>),

\begin{equation}\label{eq9}
{p_{{n^{*}}}} - {p_{{n^{'}}}} = \frac{{{G_{{n^{*}}}} - {K_{{n^{*}},{n^{'}}}}}}{L} + \sum\limits_{m = 1}^M {\left( {{c_{m,{n^{*}}}} - {c_{m,{n^{'}}}}} \right){n_m}},
\end{equation}

we can deduce that $\mathop  \cap \limits_{{n^{'}} \ne {n^{*}}}  {{p_{{n^{*}}}} > {p_{{n^{'}}}}} $ implies

\begin{equation}
{{\mathbf{\tilde A}}_{{n^{*}}}}{\mathbf{n}} \prec {{\mathbf{\tilde t}}_{{n^{*}}}},
\end{equation}

where ${{\mathbf{\tilde A}}_{{n^{*}}}}$ and ${{\mathbf{\tilde t}}_{{n^{*}}}}$ are respectively the results of removing the ${n^{*}}$-th column of ${{\mathbf{A}}_{{n^{*}}}} \in {\mathbb{R}^{N \times M}}$ and the ${n^{*}}$-th row of ${{\mathbf{t}}_{{n^{*}}}} \in {\mathbb{R}^{N \times 1}}$, with the $(n,m)$-th element of ${{\mathbf{A}}_{{n^{*}}}}$ being ${c_{m,n}} - {c_{m,{n^{*}}}}$, ${{\mathbf{t}}_{{n^{*}}}} = {\left[ {\frac{{{G_{{n^{*}}}} - {K_{{n^{*}},1}}}}{L},\frac{{{G_{{n^{*}}}} - {K_{{n^{*}},2}}}}{L}, \ldots ,\frac{{{G_{{n^{*}}}} - {K_{{n^{*}},N}}}}{L}} \right]^{\text{T}}}$, and ${\mathbf{n}} = {\left[ {{n_1},{n_2}, \ldots ,{n_M}} \right]^{\text{T}}}$.

It should be noted that although both ${{\mathbf{\tilde A}}_{{n^{*}}}}$ and ${{\mathbf{\tilde t}}_{{n^{*}}}}$ are known variables and the elements in $\bf n$ follow mutually independent Gaussian distribution, the distribution of the random variable ${{\mathbf{\tilde A}}_{{n^{*}}}}{\mathbf{n}}$ does not exist mathematically. This is because in hash beam training, it is usually the case that $N>M$, resulting in the covariance matrix ${\mathbf{\Sigma }} = {{\mathbf{\tilde A}}_{{n^{*}}}}{\mathbf{\tilde A}}_{{n^{*}}}^{\text{H}}$ not being full-rank [11]. By adding a regularization term, it can be approximated that ${{\mathbf{\tilde A}}_{{n^{*}}}}{\mathbf{n}} \sim \mathcal{N}\left( {{\mathbf{0}},{\mathbf{\Sigma }} + \lambda {\mathbf{I}}} \right)$, and the probability of successful beam alignment $\mathbb{P}\left\{ {{n^{*}} = \hat n} \right\}$ corresponds to the cumulative distribution function (CDF) value at ${{\mathbf{\tilde t}}_{{n^{*}}}}$. Unfortunately, the calculation of CDF value can only be done numerically, which prevents further analysis of the impact of codebook parameters on beam alignment performance.

§.§ Heuristic Beam Alignment Performance Metrics

To determine the optimal codebook parameter, i.e., the number of DFT beams $L$ selected by each training beam, this subsection proposes a heuristic metric for evaluating beam alignment performance. Since a fixed codebook can be optimized through extensive offline simulations, the focus here is on identifying the optimal parameter  for random codebooks. The novel metric $\tilde p$ we proposed in this study is the probability that the measurement result ${p_{{n^{*}}}}$ of the optimal DFT beam exceeds that of any other randomly selected beam measurement result ${p_{{n^{'}} \ne {n^{*}}}}$ when generating a random hash codebook.

According to (<ref>), when the hash codebook $\mathbf C$ is given,
$${p_{{n^{*}}}} - {p_{{n^{'}}}}\sim\mathcal{N}\left( {\frac{{{G_{{n^{*}}}} - {K_{{n^{*}},{n^{'}}}}}}{L},\left( {{G_{{n^{*}}}} + {G_{{n^{'}}}} - 2{K_{{n^{*}},{n^{'}}}}} \right){\sigma ^2}} \right),$$
we have

\begin{eqnarray}\label{approx problem}
\begin{aligned}
& \mathbb{P}\left\{ {{p_{{n^{*}}}} > {p_{{n^{'}}}}|{\mathbf{C}}} \right\} \\
& = \frac{1}{2}\left\{ {1 - {\text{erf}}\left( {\frac{{{G_{{n^{*}}}} - {K_{{n^{*}},{n^{'}}}}}}{{L\sqrt {2\left( {{G_{{n^{*}}}} + {G_{{n^{'}}}} - 2{K_{{n^{*}},{n^{'}}}}} \right){\sigma ^2}} }}} \right)} \right\},
\end{aligned}
\end{eqnarray}

where the values of ${G_{{n^{*}}}}$ and ${K_{{n^{*}},{n^{'}}}}$ are determined by the codebook parameter $L$. Next, we will discuss both existing  and the proposed random hash codebooks.

§.§.§ Existing random hash codebook

Utilizing combinatorial mathematics and classical probability theory, the probability of DFT beam ${n^{*}}$ being selected ${g_{{n^{*}}}}$ times is

\begin{eqnarray}\label{eq11}
\begin{aligned}
&\mathbb{P}\left\{ {{G_{{n^{*}}}} = {g_{{n^{*}}}}} \right\} = C_M^{{g_{{n^{*}}}}}{\left( {\frac{L}{N}} \right)^{{g_{{n^{*}}}}}}{\left( {1 - \frac{L}{N}} \right)^{M - {g_{{n^{*}}}}}},\\
&	{g_{{n^{*}}}} = 0,1, \dots ,M.
\end{aligned}
\end{eqnarray}

Given ${G_{{n^{*}}}} = {g_{{n^{*}}}}$, the probability of DFT beams ${n^{*}}$ and ${n^{'}}$ both occurring ${g_{{n^{*}}}}$ times is

\begin{eqnarray}
\begin{aligned}
&\mathbb{P}\left\{ {{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}}|{G_{{n^{*}}}} = {g_{{n^{*}}}}} \right\} =\\
&C_{{g_{{n^{*}}}}}^{{k_{{n^{*}},{n^{'}}}}}{\left( {\frac{{L - 1}}{{N - 1}}} \right)^{{k_{{n^{*}},{n^{'}}}}}}{\left( {1 - \frac{{L - 1}}{{N - 1}}} \right)^{{g_{{n^{*}}}} - {k_{{n^{*}},{n^{'}}}}}},\\
&	{k_{{n^{*}},{n^{'}}}} = 0,1, \ldots ,{g_{{n^{*}}}}.
\end{aligned}
\end{eqnarray}

Under the aforementioned conditions, the probability of DFT beam ${n^{'}}$ being selected ${g_{{n^{'}}}}$ times is

\begin{eqnarray}\label{eq13}
\begin{aligned}
& \mathbb{P}\left\{ {{G_{{n^{'}}}} = {g_{{n^{'}}}}|{G_{{n^{*}}}} = {g_{{n^{*}}}},{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}}} \right\} = \\
&C_{M - {k_{{n^{*}},{n^{'}}}}}^{{g_{{n^{'}}}} - {k_{{n^{*}},{n^{'}}}}}{\left( {\frac{L}{{N - 1}}} \right)^{{g_{{n^{'}}}} - {k_{{n^{*}},{n^{'}}}}}}{\left( {1 - \frac{L}{{N - 1}}} \right)^{M - {g_{{n^{'}}}}}},\\
&{g_{{n^{'}}}} = {k_{{n^{*}},{n^{'}}}},{k_{{n^{*}},{n^{'}}}} + 1, \ldots ,M.
\end{aligned}
\end{eqnarray}

Combining (<ref>)-(<ref>), we can obtain the proposed heuristic metric $\tilde{p}$ for the existing random hash codebook. The result (<ref>) and (<ref>) is shown at the top of this page.

§.§.§ Proposed random hash codebook

For the proposed random codebook, ${G_n} = G = ML/N,\forall n$, it is only necessary to analyze the distribution of  ${K_{{n^{*}},{n^{'}}}}$, which is shown in (<ref>) at the top of this page. Finally,  the proposed heuristic metric $\tilde{p}$ for the proposed random hash codebook can be derived as

\begin{equation}\label{final_result2}
\tilde p = \sum\limits_{{k_{{n^{*}},{n^{'}}}} = 0}^G {\mathbb{P}\left\{ {{p_{{n^{*}}}} > {p_{{n^{'}}}}| \mathbf{C}} \right\}\mathbb{P}\left\{ {{K_{{n^{*}},{n^{'}}}} = {k_{{n^{*}},{n^{'}}}}} \right\}}.
\end{equation}

Following the above analysis, the heuristic performance metrics  $\tilde{p}$ in (<ref>) and (<ref>) are  functions of $L$, and the optimal codebook parameters can be conveniently calculated using closed-form expressions, i.e., $L^* = \arg \mathop {\max }\limits_L \tilde p\left( L \right)$.

§ SIMULATION RESULTS

In this section, we provide numerical results to evaluate the proposed beam training method and compared with other existing schemes. We assume that the base station is equipped with $N=128$ antennas, and there are 128 DFT beams in total. The number of training beams $M=64$.

Fig. <ref> validates the theoretical analysis and shows the theoretical and simulated beam alignment performance of the proposed hash codebook and the random codebook for different $L$. We set the channel contains only LoS paths and lies in the grid point as described in Section <ref>, and the SNR was 10dB. In Fig. <ref>, the blue curve corresponds to the scale on the left side of the axis, indicating the successful probability of beam alignment $P\{n^{\ast}=\hat{n}\}$, the red curve corresponds to the logarithmic scale on the right side of the axis, and for clarity, we use the heuristic parameter as $1-\tilde{p}$. As shown in Fig. <ref>, the number of beams $L$ included in the training beam has a significant impact on the beam alignment performance. And the proposed heuristic parameter $\tilde{p}$ are approximately monotonic with the beam alignment performance. Therefore, we can obtain optimal $L^*$ by maximizing $\tilde{p}$.

Relation of $L$ to beam alignment success probability and heuristic parameter.

Fig. <ref> demonstrates the effect of $L$ on the performance of the random codebook for different SNR. It can be seen that the optimal $L$ varies with SNR. When SNR is high, a moderate value of $L$ leads to the best performance of beam alignment. However, when the SNR is low, the optimal value of $L$ is small. This is because $L$ and SNR jointly determine the mean and variance of the Gaussian random variable $p_{n^\ast}-p_{n'}$, as shown in (<ref>). In addition, when the parameter $L$ is properly set, the beam alignment performance of the proposed random codebook is much better than that of the completely random codebook.

Relation of $L$ to beam alignment success probability in different SNR.

Fig. <ref> shows beam alignment performance of different schemes in practical terahertz channels. There are $P=3$ paths in total, containing one LoS path and two NLoS paths, where the LoS path gain $\alpha_1 \sim \mathcal{CN}(0,1)$, NLoS path gain $\alpha_2,\alpha_3\sim \mathcal{CN}(0,0.01)$. The angle of the $p$-th path $\theta_p$ follows a uniform distribution in $[-\pi/2,\pi/2]$, and the fixed codebook used for simulation is selected from 1000 random codebooks. It can be seen that the beam alignment performance of the proposed hash codebook is much better than that of the completely random codebook in references [8]. At high SNR, the proposed hash codebook approaches the performance of beam sweeping, but requires only half the training overhead. In addition, compared with the proposed random hash codebook, the proposed fixed hash codebook can further improve the success rate of beam alignment by about 1%, and has the advantage of easy implementation.

Beam training success rate with $N=128, M=64, L=8, P=3$.

§ CONCLUSION

In this paper, to solve the problem of high beam training overhead for extremely large antenna arrays, we study the fast beam training scheme based on hash function and analyze the impact of hash codebook on the performance of beam alignment. We propose two codebook design methods, random and fixed, and obtain the optimal codebook generation parameters according to the theoretical analysis results. Simulation results show that the proposed scheme achieves close to the performance of beam sweeping, but only requires half of the training overhead. At the same time, the accuracy of beam alignment of the proposed scheme is much better than that of the existing schemes using completely random hash codebook.

[1]
C. J. Zhang et al., “Key technology for 5G new radio,” IEEE
Commun. Mag., vol. 56, no. 3, pp. 10–11, Mar. 2018.

[2]
L. Dai, J. Tan, Z. Chen, and H. V. Poor, “Delay-phase precoding for wideband
thz massive MIMO,” IEEE Trans. Wireless Commun., vol. 21, no. 9,
pp. 7271–7286, Sep. 2022.

[3]
J. Wang et al., “Beam codebook based beamforming protocol for
multi-gbps millimeter-wave WPAN systems,” IEEE J. Sel. Areas
Commun., vol. 27, no. 8, pp. 1390–1399, Oct. 2009.

[4]
J. Zhang, Y. Huang, Q. Shi, J. Wang, and L. Yang, “Codebook design for beam
alignment in millimeter wave communication systems,” IEEE Trans.
Commun., vol. 65, no. 11, pp. 4980–4995, Nov. 2017.

[5]
C. Qi, K. Chen, O. A. Dobre, and G. Y. Li, “Hierarchical codebook-based
multiuser beam training for millimeter wave massive MIMO,” IEEE
Trans. Wireless Commun., vol. 19, no. 12, pp. 8142–8152, Dec. 2020.

[6]
Z. Marzi, D. Ramasamy, and U. Madhow, “Compressive channel estimation and
tracking for large arrays in mm-wave picocells,” IEEE J. Sel. Topics
Signal Process., vol. 10, no. 3, pp. 514–527, Apr. 2016.

[7]
H. Yan and D. Cabric, “Compressive initial access and beamforming training for
millimeter-wave cellular systems,” IEEE J. Sel. Topics in Signal
Process., vol. 13, no. 5, pp. 1151–1166, Sep. 2019.

[8]
X. Wang, X. Hou, L. Chen, S. Suyama, and T. Asai, “Hash-based fast beam
alignment for 6G sub-terahertz MIMO,” in 2022 27th Asia Pacific
Conference on Communications (APCC), Oct. 2022, pp. 317–322.

[9]
C. You, B. Zheng, and R. Zhang, “Fast beam training for IRS-assisted
multiuser communications,” IEEE Wireless Commun. Lett., vol. 9,
no. 11, pp. 1845–1849, Nov. 2020.

[10]
Z. Li and C. Qi, “Beam training with limited feedback for multiuser mmwave
massive MIMO,” in GLOBECOM 2020 - 2020 IEEE Global Communications
Conference, Dec. 2020, pp. 1–6.

[11]
C. Raphael, “Bayesian networks with degenerate gaussian distributions,”
Methodology and Computing in Applied Probability, vol. 5, pp.
235–263, Jun. 2003.


# Quantum Annealing for Computer Vision Minimization Problems

Shahrokh Heidari  Michael J. Dinneen  Patrice Delmas School of Computer
Science, The University of Auckland, New Zealand

###### Abstract

Computer Vision (CV) labelling algorithms play a pivotal role in the domain of
low-level vision. For decades, it has been known that these problems can be
elegantly formulated as discrete energy minimization problems derived from
probabilistic graphical models (such as Markov Random Fields). Despite recent
advances in inference algorithms (such as graph-cut and message-passing
algorithms), the resulting energy minimization problems are generally viewed
as intractable. The emergence of quantum computations, which offer the
potential for faster solutions to certain problems than classical methods, has
led to an increased interest in utilizing quantum properties to overcome
intractable problems. This study investigates a new Quantum Annealing based
inference algorithm for CV discrete energy minimization problems. Our
contribution is focused on Stereo Matching as a significant CV labeling
problem. As a proof of concept, we also use a hybrid quantum-classical solver
provided by D-Wave System to compare our results with the best classical
inference algorithms in the literature.

##  1 Introduction

Computer Vision (CV) is a field of study focusing on how computers gain high-
level perception from digital images or videos, which can help decision-making
in real-world environments. While humans routinely interpret the environment,
enabling computers to perceive the real world from its representation through
images remains a largely unsolved problem. Many problems in CV are formulated
as labeling problems. A CV labeling problem consists of a set of image
features (such as pixels, edges, or image segments) on which we want to
estimate quantities from a set of labels [46] (such as intensity in Image
Restoration or disparity in Stereo Matching and Motion). Generally, CV
labeling problems are modeled by a discrete minimization problem, where an
objective function is defined to be optimized over a set of possible
labelings. When this objective function measures the badness, the optimization
problem is often called energy minimization, and the objective function is
referred to as an energy function [20]. Given the intrinsically tricky nature
of CV minimization problems, researchers have always been looking for
efficient algorithms to approximate the optimal solution as fast and
accurately as possible. Thus, there has been significant development in
minimization algorithms for CV problems from the classical methods in the 90s,
such as Simulated Annealing [12], Mean-field Annealing [21], and Iterated
Conditional Modes (ICM) [5] to the recent state-of-the-art algorithms, such as
graph-cut based [8, 6, 9, 29, 50, 51] and message-passing based [40, 19, 54]
approaches (we refer interested readers to the most recent comparative studies
on CV minimization algorithms [45, 43, 28, 44, 25]). Despite being extensively
researched and even considering the most recent advances using deep learning-
based strategies [49], which are computationally expensive, CV labeling
problems are still considered open problems with no prefect (optimal)
solutions.

Therefore, researchers have always been looking for alternatives to tackle the
problem. With the advent of quantum computations which promise potentially
lower-time complexity on certain problems than the best-classical counterparts
[14, 26, 53], recent studies have focused on leveraging quantum properties to
overcome intractable classical problems using Quantum Annealing (QA). D-Wave
System was the first company to build a Quantum Processing Unit (QPU) that
naturally approximates the ground state of a particular problem
representation, namely Ising model [38]. The importance of Ising models is
that one can solve a variety of NP-hard optimization problems by finding the
corresponding ground state [36, 10, 11]. Despite the promising experiments
[14, 26], D-Wave QPUs are specifically designed to solve optimization
problems, making them less versatile than other quantum computing approaches.
This restricts their application domain primarily to optimization and sampling
tasks, while they may not be suitable for more general-purpose computing
requirements. Also, D-Wave QPUs exhibit limited qubit connectivity, and the
scarcity of available qubits has been consistently challenging, from the
128-qubit D-Wave One built in 2011 to the newly released 5000-qubit D-Wave
Advantage. Therefore, large CV problems involving highly non-convex functions
in a search space of many thousands of dimensions have not been widely
explored to see if QA can provide advantages in real-world CV problems. In
recent years, there has been a growing interest in Quantum Computer Vision
(QCV), largely fueled by recent advancements in D-Wave QPU architectures and
their capabilities in solving optimization problems, such as Classification
[1, 15, 32, 31, 30], Synchronization [7, 2], Tracking [56], Fitting [16, 18],
Detection [34], and Matching [13, 24, 3, 55, 4] problems. However, each method
employs a distinct quantum model to represent the respective CV problem,
allowing it to be minimized on a D-Wave QPU. A versatile framework for
converting a CV problem into an appropriate quantum model holds significant
value. Such a flexible solution not only simplifies the process of adapting
various CV problems for quantum computation but also opens up new avenues for
harnessing the power of quantum computing in addressing intricate optimization
tasks.

In this study, we aim to focus on a challenging labeling problem, Stereo
Matching, and provide a general-purpose quantum model that can be used for any
CV labeling problem (such as Image Segmentation, Image Restoration, Image
Registration, Optical Flow, Object Detection, and Image Inpainting). Due to
the scarcity of available qubits on the current D-Wave QPUs, we use a D-Wave
hybrid quantum-classical solver to show the feasibility of the proposed
quantum model once enough qubits are available. Our findings show that Quantum
Annealing can offer promising results in CV applications compared to the
state-of-the-art CV minimization inference algorithms.

The paper is organized as follows: In Section 2, we provide a brief
introduction to Stereo Matching, an important CV labeling problem. In Section
3, we shift the focus to Quantum Annealing and the D-Wave QPU. Our general-
purpose quantum solution to Stereo Matching and its proof of correctness are
presented in Section 4. We provide experimental results and numerical
evaluation in Section 5. Finally, Section 6 concludes the paper.

##  2 Stereo Matching

The characteristics of binocular vision in humans allow for the simultaneous
observation of a singular object by both eyes. This ability significantly
contributes to the understanding of depth in the brain. The distance between
our eyes, often referred to as “baseline”, facilitates slight variation in the
perspective captured by each eye. Despite each eye observing a nearly
identical image, a marginal displacement exists. The brain uses this
displacement to perceive a 3D observation from the scene. Likewise, a stereo
vision system is designed to replicate human vision mechanisms. This system
comprises two horizontal cameras on the left and right sides, effectively
simulating human binocular perception. Each camera in the system records an
image that, while fundamentally similar, features a certain degree of
displacement. This displacement, often called disparity, signifies the
difference in the position of a 3D point, as observed from two different
viewpoints (the left and right viewpoints) [22]. By implementing a stereo
vision system, the main goal is to construct a 3D model using the left and
right stereo images. This procedure may encompass various stages, including
Camera Calibration (optional), Rectification, Stereo Matching, and 3D
Reconstruction [37]. Identifying corresponding 2D projections in stereo images
is known as Stereo Matching. This process represents the most computationally
demanding component of a stereo vision system. When the left image serves as
the reference, the goal of Stereo Matching is to find the corresponding pixel
in the right image for every pixel in the left image. Each matched pixel pair
represent projections of a 3D point in the real world. Stereo Matching methods
are broadly categorized into global and local approaches. While local methods
prioritize speed, often at the cost of accuracy due to susceptibilities like
local ambiguities and occlusions, global methods comprehensively consider the
entire image during disparity computation. Although computationally demanding,
they effectively address challenges such as occluded and textureless regions
[22]. These methods typically lean on probabilistic graphical models, a potent
blend of probability and graph theory, for their formalism [52]. Based on the
defined probabilistic graphical model, an energy function is modeled which can
be minimized to solve the Stereo Matching problem [52, p. 1612]. In the
following, we provide the general form of a global Stereo Matching energy
function, which can be adapted for any CV labeling problem (see the recent
comparative study on CV labeling problems [25] for more information).

Let $I_{l}$ and $I_{r}$ be a pair of $n\times m$ stereo images, and
$D=\\{d_{min},\dots,d_{max}\\}$ be a set of positive integers, where $d_{min}$
and $d_{max}$ are the lowest and highest possible disparity values,
respectively. Considering the left image $I_{l}$ as the reference for which we
want to compute a disparity map, the set of pixels is defined as (1). Here, we
also initialize $N$ as a 4-neighborhood system defined in (2).

$\small P=\\{(i,j)\mid i\in\\{0,\dots,n-1\\},j\in\\{0,\dots,m-1\\}\\},$ (1)
$\displaystyle N$
$\displaystyle=\\{\\{(i,j),(i^{\prime},j^{\prime})\\}\mid(i,j)\in P,$ (2)
$\displaystyle(i^{\prime},j^{\prime})\in\\{(i\pm 1,j),(i,j\pm 1)\\}\\},$

where, $0\leq i^{\prime}<n,0\leq j^{\prime}<m$. In a global Stereo Matching
model, the Stereo Matching problem is modeled by a labeling problem where each
pixel in $P$ is labeled by a disparity value in $D$ [46, p. 5]. In fact, a
“labeling” involves mapping from $P$ to $D$. Such a labeling problem is
defined by a discrete optimization problem, where an energy function is
defined to be minimized over a set of possible labelings. This energy function
has two terms. The first term penalizes the solutions when they are
inconsistent with the data, and the second term imposes some constraints on
spatial coherence [44, p. 1]. Let $\mathbf{w}\in D^{n\times m}$ be a vector of
variables defined as $\mathbf{w}=(w_{i,j})_{(i,j)\in P}$, where $w_{i,j}\in
D$. The global Stereo Matching energy function $F:D^{n\times
m}\rightarrow\mathbb{R}^{+}$ is defined as (3).

$\displaystyle F(\mathbf{w})$ $\displaystyle=\sum_{(i,j)\in
P}\theta_{\\{i,j\\}}(w_{i,j})+\lambda\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in
N}\delta(w_{i,j},w_{i^{\prime},j^{\prime}})$ (3)

where,

$\displaystyle\theta_{i,j}(w_{i,j})=\left|I_{l}(i,j)-I_{r}(i-w_{i,j},j)\right|,$
$\displaystyle\delta(w_{i,j},w_{i^{\prime},j^{\prime}})=\left\\{\begin{array}[]{ll}0,&\hbox{if
$w_{i,j}=w_{i^{\prime},j^{\prime}}$;}\\\ 1,&\hbox{otherwise.}\\\
\end{array}\right.$

The first term is the Sum of Absolute Difference (SAD) matching cost function
defined by $\theta_{\\{i,j\\}}:D\rightarrow\mathbb{R}^{+}$. When
$\theta_{\\{i,j\\}}(w_{i,j})$ is (or close to) zero, it means the pixel
$(i,j)$ in the left image matches the pixel $(i-w_{i,j},j)$ in the right
image, and they are more likely to be the same projections of a 3D point in
the real world. In the second term, $\delta:D^{2}\rightarrow\\{0,1\\}$ is the
penalty function that penalizes the variation of the disparities, adding one
when the allocated disparities to a pair of neighboring pixels are not equal
and zero otherwise. The second term assumes that the disparities of a
neighborhood of pixels present some coherence and generally do not change
abruptly [35]. Furthermore, $\lambda\in\mathbb{R}^{+}$, known as the
smoothness factor, weighs the penalties given by the second term.

We aim to provide a general-purpose quantum model for the defined global
Stereo Matching problem, which can be used for any CV labeling problem. Thus,
we first give the preliminaries to describe this quantum model.

##  3 Quantum Annealing

Quantum Annealing [17] is a specialized optimization technique that leverages
principles from quantum mechanics to solve complex computational problems. In
this model, quantum bits (qubits) are particles in a quantum dynamical system
that evolve based on special forces acting on them. These forces are either
internal (from interactions among qubits) or external (from other sources).
Each state of a register of qubits has energy based on the applied forces. A
time-dependent Hamiltonian is a mathematical representation of a system,
providing information about the system’s energy and detailing the forces
acting upon it at any given time [38]. Quantum Annealing is a computational
technique employed to discover the state of the system with the minimum energy
as determined by the time-dependent Hamiltonian. Consequently, Quantum
Annealing constitutes a computational paradigm known for its efficiency in
addressing optimization problems and providing approximations to the optimal
solutions. It is inspired by the concept of annealing in metallurgy, where a
material is slowly cooled to minimize defects and reach a low-energy state. In
Quantum Annealing, this cooling process is simulated by a QPU known as a
quantum annealer which is based on a time-dependent Hamiltonian $H(t)$ that
has three components [38]: Initial Hamiltonian $H_{I}$, where all qubits are
in a superposition state. Problem Hamiltonian $H_{p}$, where the specific
forces are defined to encode the objective function. The lowest-energy state
of $H_{p}$ is the solution that minimizes the objective function. Adiabatic
path $s(t)$, which is a smooth function that decreases from $1$ to $0$, such
as $s(t)=1-\frac{t}{t_{f}}$, where $s(t)$ decreases from $1$ to $0$ as $t$
increases from $0$ to some elapsed time $t_{f}$. During Quantum Annealing, the
Initial Hamiltonian is slowly evolved along the Adiabatic path to the Problem
Hamiltonian as $H(t)=s(t)H_{I}+(1-s(t))H_{p}$ [38], decreasing the influence
of $H_{I}$ over time to reach $H_{P}$ as $s(t)$ goes from $1$ to $0$. D-Wave
Systems was the first company to build a quantum annealer. To
minimize/maximize an objective function using Quantum Annealing and a D-Wave
QPU, it should be in a standard model like Ising or Quadratic Unconstrained
Binary Optimization (QUBO) models [38]. Given a vector of $n$ binary variables
as $\mathbf{x}=(x_{1},x_{2},\dots,x_{n})\in\\{0,1\\}^{n}$, a QUBO model is
represented as $H_{qubo}(\mathbf{x})=\mathbf{x}^{T}\mathbf{Q}\mathbf{x}$,
where $\\{0,1\\}^{n}$ is a set of $n$ binary values, and $\mathbf{Q}$ is an
$n\times n$ matrix that can be chosen to be upper-diagonal. Therefore,
$H_{qubo}(\mathbf{x})$ can be reformulated as (4).

$H_{qubo}(\mathbf{x})=\sum_{i}\mathbf{Q}_{i,i}x_{i}+\sum_{i<j}\mathbf{Q}_{i,j}x_{i}x_{j}.$
(4)

The diagonal terms $\mathbf{Q}_{i,i}$ are the linear coefficients acting as
the external forces, and the off-diagonal terms $\mathbf{Q}_{i,j}$ are the
quadratic coefficients for the internal forces [38].

##  4 Quantum Stereo Matching

We introduce an equivalent QUBO model to the global Stereo Matching
minimization problem (3) and provide proof of its validity. Our idea draws
inspiration from the approach employed by the D-Wave Ocean SDK when handling
discrete objective functions [42]. We first allocate $|D|$ binary variables to
each pixel $(i,j)\in P$, where $|D|$ is the number of elements in $D$, $0\leq
i\leq n-1$, and $0\leq j\leq m-1$. Therefore, we define
$\mathbf{x}\in\\{0,1\\}^{nm|D|}$ as a vector of $nm|D|$ binary variables such
that $\mathbf{x}=(x_{i,j,d})$ for all $(i,j)\in P$ and $d\in D$. Let our QUBO
model be defined as (5).

$\displaystyle H(\mathbf{x})$ $\displaystyle=\alpha\sum_{(i,j)\in
P}\left(1-\sum_{d\in D}x_{i,j,d}\right)^{2}+\sum_{(i,j)\in P}\sum_{d\in
D}\theta_{\\{i,j\\}}(d)x_{i,j,d}$ (5)
$\displaystyle+\lambda\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in
N}\sum_{d_{1}\in D}\sum_{d_{2}\in
D}\delta(d_{1},d_{2})x_{i,j,d_{1}}x_{i^{\prime},j^{\prime},d_{2}},$

where $\alpha>\left(\sum_{(i,j)\in P}\max\\{\theta_{\\{i,j\\}}(d)\mid d\in
D\\}\right)+\lambda|N|$, and $|N|$ is the number of elements in $N$. We set
$\mathbf{x}^{*}=\operatorname*{arg\,min}_{\mathbf{x}}H(\mathbf{x})$ and define
a vector of $nm$ integer values as $\mathbf{w}^{*}=(w^{*}_{i,j})_{(i,j)\in
P}$, where $w^{*}_{i,j}=d$ if $x^{*}_{i,j,d}=1$. Then, $\mathbf{w}^{*}$
minimizes the global Stereo Matching energy function (3).

### Proof of Correctness

Eq. 5 has three parts. The first part guarantees each pixel is assigned a
unique disparity value from $D$. The second calculates the cost of the
assigned disparity values to the pixels. The third part encodes the defined
contextual constraint.

######  Definition 1.

$\mathbf{x}$ is called feasible if and only if $\sum_{d\in D}x_{i,j,d}=1$ for
all pixels $(i,j)\in P$. We denote a feasible $\mathbf{x}$ by
$\mathbf{x}^{\prime}$.

Definition 1 states that given a pixel $(i,j)\in P$, its corresponding vector
of binary variables
$(x^{\prime}_{i,j,d_{min}},\dots,x^{\prime}_{i,j,d_{max}})$ has only one value
of “1” in its values, making it possible to label each pixel uniquely by a
disparity $d\in D$. Hence, the allocated disparity to a pixel $(i,j)\in P$ is
$d$ if $x^{\prime}_{i,j,d}=1$.

######  Definition 2.

Given $\mathbf{x}^{\prime}$, the corresponding integer vector
$\mathbf{w}^{\prime}=(w^{\prime}_{i,j})_{(i,j)\in P}$ is called a _labeling_ ,
where $w^{\prime}_{i,j}=d$ if $x^{\prime}_{i,j,d}=1$.

######  Lemma 1.

Given a feasible $\mathbf{x}^{\prime}$ and its corresponding labeling
$\mathbf{w}^{\prime}$, the equality
$H(\mathbf{x}^{\prime})=F(\mathbf{w}^{\prime})$ holds, where $F$ is the global
Stereo Matching energy function in (3).

###### Proof.

Considering $H(\mathbf{x}^{\prime})$ in (5),

  * •

Since $\mathbf{x}^{\prime}$ is feasible, $\sum_{d\in D}x^{\prime}_{i,j,d}=1$
for all pixels $(i,j)\in P$ by Definition 1. Therefore, we have

$\alpha\sum_{(i,j)\in P}\left(1-\sum_{d\in
D}x^{\prime}_{i,j,d}\right)^{2}=\alpha\sum_{(i,j)\in
P}\left(1-1\right)^{2}=0.$

  * •

Given a pixel $(i,j)\in P$, only one variable in the vector
$(x^{\prime}_{i,j,d_{min}},\dots,x^{\prime}_{i,j,d_{max}})$ is one, and all
the others are zero. This non-zero variable is
$x^{\prime}_{i,j,w^{\prime}_{i,j}}$ by Definition 2. Therefore, we have

$\displaystyle\sum_{(i,j)\in P}\sum_{d\in
D}\theta_{\\{i,j\\}}(d)x^{\prime}_{i,j,d}$ $\displaystyle=\sum_{(i,j)\in
P}\theta_{\\{i,j\\}}(w^{\prime}_{i,j})x^{\prime}_{i,j,w^{\prime}_{i,j}}=\sum_{(i,j)\in
P}\theta_{\\{i,j\\}}(w^{\prime}_{i,j}).$

  * •

Given $\\{(i,j),(i^{\prime},j^{\prime})\\}\in N$, the two corresponding
vectors of binary variables are

    * –

$(i,j):(x_{i,j,d_{min}},\dots,x_{i,j,d_{max}})$

    * –

$(i^{\prime},j^{\prime}):(x^{\prime}_{i^{\prime},j^{\prime},d_{min}},\dots,x^{\prime}_{i^{\prime},j^{\prime},d_{max}})$

Since $\mathbf{x}^{\prime}$ is feasible, only one of the variables in each
vector is one, and the others are zero. These variables are
$x^{\prime}_{i,j,w^{\prime}_{i,j}}$ and
$x^{\prime}_{i^{\prime},j^{\prime},w^{\prime}_{i^{\prime}j^{\prime}}}$,
respectively, by Definition 2. Thus, we can write

$\displaystyle\lambda\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in
N}\sum_{d_{1}\in D}\sum_{d_{2}\in
D}\delta(d_{1},d_{2})x^{\prime}_{i,j,d_{1}}x^{\prime}_{i^{\prime},j^{\prime},d_{2}}$
$\displaystyle=\lambda\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in
N}\delta(w^{\prime}_{i,j},w^{\prime}_{i^{\prime}j^{\prime}})x^{\prime}_{i,j,w^{\prime}_{i,j}}x^{\prime}_{i^{\prime},j^{\prime},w^{\prime}_{i^{\prime}j^{\prime}}}$
$\displaystyle=\lambda\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in
N}\delta(w^{\prime}_{i,j},w^{\prime}_{i^{\prime}j^{\prime}}).$

Therefore, we can rewrite $H(\mathbf{x}^{\prime})$ as follows.

$\displaystyle H(\mathbf{x}^{\prime})$ $\displaystyle=\sum_{(i,j)\in
P}\theta_{\\{i,j\\}}(w^{\prime}_{i,j})+\lambda\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in
N}\delta(w^{\prime}_{i,j},w^{\prime}_{i^{\prime}j^{\prime}})=F(\mathbf{w}^{\prime}).$

∎

######  Lemma 2.

Let $\mathbf{x}^{*}=\operatorname*{arg\,min}_{\mathbf{x}}H(\mathbf{x})$.
$\mathbf{x}^{*}$ is feasible.

###### Proof.

For ease of reference, we rewrite $H(\mathbf{x})$ as follows:

$H(\mathbf{x})=\alpha\mathcal{A}(\mathbf{x})+\mathcal{B}(\mathbf{x}),$

where

$\displaystyle\mathcal{A}(\mathbf{x})$ $\displaystyle=\sum_{(i,j)\in
P}\left(1-\sum_{d\in D}x_{i,j,d}\right)^{2}$
$\displaystyle\mathcal{B}(\mathbf{x})$ $\displaystyle=\sum_{(i,j)\in
P}\sum_{d\in
D}\theta_{\\{i,j\\}}(d)x_{i,j,d}+\lambda\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in
N}\sum_{d_{1}\in D}\sum_{d_{2}\in
D}\delta(d_{1},d_{2})x_{i,j,d_{1}}x_{i^{\prime},j^{\prime},d_{2}}.$

Towards a contradiction, suppose that $\mathbf{x}^{*}$ is not feasible. In
this case, $\mathcal{A}(\mathbf{x}^{*})\neq 0$ and it is non-negative.
Therefore,

$H(\mathbf{x}^{*})=\alpha\mathcal{A}(\mathbf{x}^{*})+\mathcal{B}(\mathbf{x}^{*}).$
(6)

Given a feasible $\mathbf{x}^{\prime}$, $\mathcal{A}(\mathbf{x}^{\prime})=0$,
and we have

$H(\mathbf{x}^{\prime})=\mathcal{B}(\mathbf{x}^{\prime}).$ (7)

Since $\mathbf{x}^{\prime}$ is feasible, $\mathcal{B}(\mathbf{x}^{\prime})$
adds penalty values up to a maximum of

$\left(\sum_{(i,j)\in P}\max\\{\theta_{i,j}(d)\mid d\in
D\\}\right)+\lambda|N|.$

Considering (7), we have

$H(\mathbf{x}^{\prime})\leq\left(\sum_{(i,j)\in P}\max\\{\theta_{i,j}(d)\mid
d\in D\\}\right)+\lambda|N|.$ (8)

We know that $\alpha>\left(\sum_{(i,j)\in P}\max\\{\theta_{i,j}(d)\mid d\in
D\\}\right)+\lambda|N|$, $\mathcal{A}(\mathbf{x}^{*})$ is non-zero and non-
negative, and $\mathcal{B}(\mathbf{x}^{*})$ is non-negative. Considering (6),
we can write

$\displaystyle H(\mathbf{x}^{*})$
$\displaystyle=\alpha\mathcal{A}(\mathbf{x}^{*})+\mathcal{B}(\mathbf{x}^{*})>\left(\sum_{(i,j)\in
P}\max\\{\theta_{i,j}(d)\mid d\in D\\}\right)+\lambda|N|.$ (9)

The following statement is true by (8) and (9):
$H(\mathbf{x}^{\prime})<H(\mathbf{x}*)$, which is a contradiction. Therefore,
$\mathbf{x}^{*}$ is feasible. ∎

######  Theorem 1.

Given $\mathbf{w}^{*}$ as the corresponding labeling of $\mathbf{x}^{*}$,
$\mathbf{w}^{*}$ minimizes the global Stereo Matching energy function $F$
defined in (3).

###### Proof.

Towards a contradiction, we suppose that $\mathbf{w}^{*}$ does not minimize
$F$. In this case, there must be a feasible $\mathbf{x^{\prime}}$ for which
its corresponding labeling $\mathbf{w}^{\prime}$ minimizes $F$. Therefore, we
have $F(\mathbf{w}^{\prime})<F(\mathbf{w}^{*})$. Since $\mathbf{x}^{\prime}$
and $\mathbf{x}^{*}$ are both feasible (see Lemma 2), we have
$H(\mathbf{x}^{\prime})<H(\mathbf{x}^{*})$ by Lemma 1. This is a contradiction
because in this case
$\mathbf{x}^{*}\neq\operatorname*{arg\,min}_{\mathbf{x}}H(\mathbf{x})$. ∎

Eq. 5 is versatile and can be adapted for a variety of CV labeling problems by
replacing $P$ with any desired set of image features, and replacing $D$ with
an appropriate set of labels depending on the application. Then, the first and
second terms in (3) can be defined accordingly. The modified QUBO remains
consistent with the QUBO model described in (5). This adaptability showcases
the broader applicability of the model, making it a flexible tool for
addressing a range of CV labeling challenges.

##  5 Evaluation and Experimental Results

###  5.1 Qubit Complexity

D-Wave quantum computers have showcased remarkable potential in solving
optimization problems. However, one significant challenge they face is the
limited availability of qubits. D-Wave QPUs employ Quantum Annealing to find
(estimate) the global minimum of a QUBO/Ising model. While effective for
specific problem types, this approach often requires a large number of qubits,
and the current generation of D-Wave QPUs have constraints on the number of
qubits that can be utilized. Consequently, proposing a QUBO model with fewer
variables is paramount as it addresses the current limitations in qubit
availability, enables the solution of larger and more complex problems, widens
access to Quantum Annealing, and enhances the robustness and practicality of
Quantum Annealing technology in solving real-world optimization challenges.
Recall $P$ as the set of pixels for a pair of stereo images with size $n\times
m$, and $D$ as the set of possible disparities values, where $|P|=nm$ and
$|D|=k$ denotes the number of elements in $P$ and $D$, respectively. Given the
defined vector of binary variables in (5), the number of QUBO variables in our
general-purpose quantum model is $nmk$. Table 1 compares our quantum model
with the existing labeling-based quantum solutions that can be utilized for
Stereo Matching. Cruz-Santos et al. [13] and Heidari et al.[24] models are
based on the minimum st cut problem, and Heidari et al.[23] approach reduces a
CV labeling problem to the minimum multi-way cut problem.

Table 1: Qubit-complexity comparison of the proposed quantum Stereo Matching models. Model | Qubit complexity
---|---
[13] | $7nmk+9nm-2nk-2mk-2n-2m+2$
[24] | $nmk+nm+2$
[23] | $nmk+k^{2}$
Ours | $nmk$

###  5.2 Experimental Results

Once a QUBO model is prepared, it needs to be embedded within the QPU hardware
architecture for the minimization process. Embedding is the crucial step of
mapping QUBO variables onto the available qubits on the hardware. Embedding
can be challenging due to the relatively limited qubits and the restricted
hardware connectivity. Consequently, it is common to chain two or more qubits
together on the QPU to represent a single QUBO variable. While many real-world
applications can successfully run on the D-Wave QPUs, there are cases where
the input data is too large to be directly solved by Quantum Annealing,
primarily because of the qubit scarcity. To overcome this size limitation,
hybrid solvers combine classical and quantum approaches for problem-solving.
D-Wave hybrid solvers can handle problems with a significantly higher number
of variables than those directly solvable by a D-Wave QPU, offering a reliable
estimate of the future accuracy of D-Wave QPUs once more qubits become
available on the hardware. As a proof of concept, we utilize the Constrained
Quadratic Model (CQM) D-Wave hybrid solver to minimize the proposed quantum
Stereo Matching model. This solver has the capability of handling up to
500,000 QUBO variables, but it still poses restrictions on the size of the
input stereo images and the number of disparities that can be processed.
Therefore, we had to use cropped pairs of stereo images to analyze the
performance of the Stereo Matching quantum model. We chose four pairs of
stereo images from 2001-Middlebury image datasets [47], namely Venus and Bull,
Sawtooth, and Barn. We could not use the latest stereo datasets because of
their high disparity range. Given a pair of cropped regions from recent
Middlebury stereo datasets, the majority of regions in both cropped stereo
images would be occluded due to a large disparity range, resulting in only a
small portion of the scene being visible in both images. This makes them not
suitable to evaluate our quantum model due to the simplicity of the defined
global Stereo Matching energy function in (3). To identify a more common
region of interest in both stereo images, we selected our pairs of stereo
images from the 2001 Middlebury image dataset [47], as well as two “natural”
images (Tree and Castle) from the real world with low disparity ranges to
incorporate complex scene structures into our prepared dataset. Figure 1
illustrates our prepared stereo dataset with the corresponding “ground
truths.” We did not have ground truths for the natural images since they were
not created in controlled laboratory settings like 2001-Middlebury image
datasets [47]. Therefore, we incorporated a deep-learning-based model [33] to
get fairly accurate disparity maps to be used as the corresponding ground
truths. Note that we used the gray-scale versions of the shown stereo images.

We also selected the best-performing and state-of-the-art classical
minimization algorithms commonly used in CV, so that we can compare our
quantum model with the classical counterparts. The selected classical
algorithms includes Swap move [9], Expansion move [9], Max product Loopy
Belief Propagation (LBP) [19], and the improved Tree Re-weighted Message
Passing (TRW-S) [27]. We utilized the Middlebury software framework for our
classical implementations [48].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: The prepared stereo datasets, (a) Venus, (b) Bull, (c) Sawtooth, (d)
Barn, (c) Tree, and (d) Castle. In each row of images, we have the left stereo
image, the right stereo image, and the corresponding ground truth for the left
stereo image. The white squares show the cropped regions for our experiments.

Considering the global Stereo Matching energy function defined in (3), we
established an initial $\lambda=20$ for all the benchmark minimization
algorithms to ensure a fair comparison. Figure 2 shows the computed disparity
maps by the benchmark minimization algorithms.

Figure 2: Computed disparity maps by the benchmark minimization algorithms.

Next, we define two widely-used metrics [39] to evaluate the accuracy of our
results in comparison to the corresponding ground truths: the root-mean-
squared error ($rms$) and the percentage of mismatched pixels (bad-$\beta$).
Given a disparity map $\mathcal{D}$ and a ground truth $\mathcal{T}$ defined
by $n\times m$ matrices of integers, $rms$ and bad-$\beta$ are defined as
(10).

$\hbox{{rms}}=\sqrt{\frac{1}{nm}\sum_{i=1}^{n}\sum_{j=1}^{m}\left(\mathcal{D}_{i,j}-\mathcal{T}_{i,j}\right)^{2}},$
(10)
$\hbox{{bad-$\beta$}}=\left(\frac{1}{nm}\sum_{i=1}^{n}\sum_{j=1}^{m}\left(|\mathcal{D}_{i,j}-\mathcal{T}_{i,j}|>\beta\right)\right)\times
100$ (11)

where $\beta\in\mathbb{R}^{+}$ is the disparity error tolerance. In the
following evaluation, we set $\beta$ to $0.5$ and $1.0$, named bad-$0.5$ an
bad-$1.0$, respectively. Table. 2 compares the performance of each
minimization algorithm based on the defined metrics. Given the cropped stereo
images, the results suggest that the quantum model outperformed the classical
counterparts on the Bull, Sawtooth. and Tree datasets, and performed
competitively on the Venus, Barn, and Castle datasets. Our findings show that
Quantum Annealing can offer promising results in CV applications compared to
the state-of-the-art CV minimization inference algorithms. Due to the scarcity
of available qubits on the current D-Wave QPUs, we were not able to use a pure
Quantum Annealing minimization, and we used a D-Wave hybrid solver, which
offers a reliable estimate of the future accuracy of D-Wave QPUs once more
qubits become available on the hardware.

Table 2: Numerical evaluation for the Venus stereo dataset. Dataset | Method | rms | bad-0.5 (%) | bad-1.0 (%)
---|---|---|---|---
Venus | Ours | 2.25 | 40.44 | 10.45
Swap | 2.09 | 47.57 | 10.23
Expansion | 1.94 | 43.81 | 9.76
BP-M | 1.96 | 47.09 | 9.30
TRW-S | 1.92 | 44.07 | 9.54
Bull | Ours | 2.33 | 36.07 | 7.08
Swap | 2.38 | 37.43 | 7.40
Expansion | 2.38 | 37.28 | 7.36
BP-M | 2.39 | 37.42 | 7.25
TRW-S | 2.38 | 37.06 | 7.32
Sawtooth | Ours | 2.27 | 22.54 | 10.26
Swap | 2.44 | 22.76 | 10.27
Expansion | 2.41 | 22.85 | 10.44
BP-M | 2.41 | 23.59 | 10.30
TRW-S | 2.36 | 22.67 | 10.36
Barn | Ours | 2.27 | 14.37 | 7.41
Swap | 2.23 | 16.11 | 7.51
Expansion | 2.21 | 16.09 | 7.55
BP-M | 2.38 | 20.21 | 8.25
TRW-S | 2.23 | 15.54 | 7.33
Tree | Ours | 2.99 | 24.99 | 13.22
Swap | 3.32 | 33.85 | 14.39
Expansion | 3.27 | 31.73 | 13.23
BP-M | 3.16 | 33.45 | 13.91
TRW-S | 3.12 | 31.81 | 13.22
Castle | Ours | 2.74 | 34.52 | 17.62
Swap | 2.83 | 32.85 | 17.25
Expansion | 2.76 | 33.68 | 16.94
BP-M | 2.99 | 41.12 | 21.21
TRW-S | 2.66 | 33.36 | 16.62

Since our model (5) is a direct equivalent to the global Stereo Matching
energy function (3), its energy solution can be compared with that of the
iterative classical minimization algorithms. Figure 3 shows the energies of
the solutions obtained by each minimization model over the provided stereo
datasets. According to the findings, our approach demonstrated a capacity to
obtain solutions of lower energy in comparison to the iterative classical
minimization methods for each provided stereo dataset. This observation
underscores the effectiveness of QUBOs when solved by D-Wave hybrid solvers.
We do not provide a comparison in terms of running time, as the classical
iterative minimization algorithms were significantly faster than the D-Wave
hybrid solver when minimizing the corresponding QUBO models. The reason is
because of the way that a D-Wave hybrid solver works. A D-Wave hybrid solver
is based on the D-Wave Hybrid Solver Service (HSS). Once a QUBO is provided to
the HSS, it activates one or more heuristic solvers that run in parallel,
either on a CPU or a GPU platform, to identify high-quality solutions. Each
heuristic solver comprises a classical heuristic module that navigates the
search space, and a quantum module is responsible for formulating quantum
queries directed to the D-Wave Advantage QPU. Solutions retrieved from the QPU
assist the heuristic modules in pinpointing more viable search space regions
or refining the current solutions. Each heuristic solver forwards its top
solution to the HSS solver. The HSS solver then determines the best solution
from the collective outputs of the heuristic solvers and relays this optimal
solution back to the user [41]. Therefore, the running time is not derived
from a direct QPU minimization to be compared with the classical minimization
methods. We used a D-Wave hybrid solver as proof of concept to show that
Quantum Annealing is capable of being used in CV labeling problems once enough
qubits are available on the QPU in future.

Figure 3: A comparison between energies obtained by our model and the
benchmarking classical minimization.

##  6 Conclusion

CV labelling algorithms play a pivotal role in the domain of low-level vision.
For decades, it has been known that these problems can be elegantly formulated
as discrete energy-minimization problems derived from probabilistic graphical
models. Despite recent advances in inference algorithms, the resulting energy-
minimization problems are generally viewed as intractable. In this study, we
presented a Quantum Annealing-based method for solving CV discrete
optimization problems, specifically for Stereo Matching. However, our proposed
quantum model is not limited to Stereo Matching and can be applied to various
CV labeling problems such as Image Segmentation, Image Restoration, Image
Registration, Optical Flow, Object Detection, and Image Inpainting. We
provided proof of correctness to demonstrate the equivalence of the proposed
quantum model to the original discrete minimization energy function. Due to
the limited availability of qubits on the quantum hardware, we were not able
to minimize the Stereo Matching energy function directly on the QPU. Instead,
we utilized a D-Wave hybrid solver to show the feasibility of our proposed
quantum model. Our results showed promising solutions with lower energies
compared to the best classical minimization algorithms in the literature. When
there are enough qubits available, it may be a subject for future research to
determine if a quantum-based CV inference offers any advantages over classical
minimization methods in terms of accuracy and speed.

## Acknowledgements

We thank Prof. Cristian Calude for his helpful discussions on this paper.

## References

  * [1] Steven H Adachi and Maxwell P Henderson.  Application of quantum annealing to training of deep neural networks.  arXiv preprint arXiv:1510.06356, 2015.
  * [2] Federica Arrigoni, Willi Menapace, Marcel Seelbach Benkner, Elisa Ricci, and Vladislav Golyanik.  Quantum motion segmentation.  In European Conference on Computer Vision, pages 506–523. Springer, 2022.
  * [3] Marcel Seelbach Benkner, Vladislav Golyanik, Christian Theobalt, and Michael Moeller.  Adiabatic quantum graph matching with permutation matrix constraints.  In 2020 International Conference on 3D Vision (3DV), pages 583–592. IEEE, 2020.
  * [4] Marcel Seelbach Benkner, Maximilian Krahn, Edith Tretschk, Zorah Lähner, Michael Moeller, and Vladislav Golyanik.  Quant: Quantum annealing with learnt couplings.  arXiv preprint arXiv:2210.08114, 2022.
  * [5] Julian Besag.  On the statistical analysis of dirty pictures.  Journal of the Royal Statistical Society: Series B (Methodological), 48(3):259–279, 1986.
  * [6] Stan Birchfield and Carlo Tomasi.  Multiway cut for stereo and motion with slanted surfaces.  In Proceedings of the seventh IEEE international conference on computer vision, volume 1, pages 489–495. IEEE, 1999.
  * [7] Tolga Birdal, Vladislav Golyanik, Christian Theobalt, and Leonidas J Guibas.  Quantum permutation synchronization.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13122–13133, 2021.
  * [8] Yuri Boykov, Olga Veksler, and Ramin Zabih.  Markov random fields with efficient approximations.  In Proceedings. 1998 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (Cat. No. 98CB36231), pages 648–655. IEEE, 1998.
  * [9] Yuri Boykov, Olga Veksler, and Ramin Zabih.  Fast approximate energy minimization via graph cuts.  IEEE Transactions on pattern analysis and machine intelligence, 23(11):1222–1239, 2001.
  * [10] Cristian S. Calude and Michael J. Dinneen.  Solving the broadcast time problem using a D-Wave quantum computer.  In Advances in Unconventional Computing, pages 439–453. Springer, 2017.
  * [11] Cristian S. Calude, Michael J. Dinneen, and Richard Hua.  QUBO formulations for the graph isomorphism problem and related problems.  Theoretical Computer Science, 701:54–69, 2017.
  * [12] Vladimír Černỳ.  Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm.  Journal of optimization theory and applications, 45(1):41–51, 1985\.
  * [13] William Cruz-Santos, Salvador E Venegas-Andraca, and Marco Lanzagorta.  A QUBO formulation of the stereo matching problem for D-Wave quantum annealers.  Entropy, 20(10):786, 2018.
  * [14] Vasil S Denchev, Sergio Boixo, Sergei V Isakov, Nan Ding, Ryan Babbush, Vadim Smelyanskiy, John Martinis, and Hartmut Neven.  What is the computational value of finite-range tunneling?  Physical Review X, 6(3):031015, 2016.
  * [15] Vivek Dixit, Raja Selvarajan, Muhammad A. Alam, Travis S. Humble, and Sabre Kais.  Training and Classification using a Restricted Boltzmann Machine on the D-Wave 2000Q.  arXiv:2005.03247 [cs, stat], May 2020.  arXiv: 2005.03247.
  * [16] Anh-Dzung Doan, Michele Sasdelli, David Suter, and Tat-Jun Chin.  A hybrid quantum-classical algorithm for robust fitting.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 417–427, 2022.
  * [17] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser.  Quantum computation by adiabatic evolution.  arXiv preprint quant-ph/0001106, 2000.
  * [18] Matteo Farina, Luca Magri, Willi Menapace, Elisa Ricci, Vladislav Golyanik, and Federica Arrigoni.  Quantum multi-model fitting.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13640–13649, 2023.
  * [19] Pedro F Felzenszwalb and Daniel P Huttenlocher.  Efficient belief propagation for early vision.  International journal of computer vision, 70(1):41–54, 2006.
  * [20] Pedro F Felzenszwalb and Ramin Zabih.  Dynamic programming and graph algorithms in computer vision.  IEEE transactions on pattern analysis and machine intelligence, 33(4):721–740, 2010.
  * [21] Davi Geiger and Federico Girosi.  Parallel and deterministic algorithms from mrfs: Surface reconstruction.  IEEE Transactions on Pattern Analysis & Machine Intelligence, 13(05):401–412, 1991.
  * [22] Rostam Affendi Hamzah and Haidi Ibrahim.  Literature survey on stereo vision disparity map algorithms.  Journal of Sensors, 2016, 2016.
  * [23] Shahrokh Heidari, Michael J Dinneen, and Patrice Delmas.  An equivalent QUBO Model to the minimum multi-way cut problem.  Technical report, Department of Computer Science, The University of Auckland, New Zealand, 2022.
  * [24] Shahrokh Heidari, Mitchell Rogers, and Patrice Delmas.  An improved quantum solution for the stereo matching problem.  In 2021 36th International Conference on Image and Vision Computing New Zealand (IVCNZ), pages 1–6. IEEE, 2021.
  * [25] Jörg H Kappes, Bjoern Andres, Fred A Hamprecht, Christoph Schnörr, Sebastian Nowozin, Dhruv Batra, Sungwoong Kim, Bernhard X Kausler, Thorben Kröger, Jan Lellmann, et al.  A comparative study of modern inference techniques for structured discrete energy minimization problems.  International Journal of Computer Vision, 115(2):155–184, 2015\.
  * [26] James King, Sheir Yarkoni, Jack Raymond, Isil Ozfidan, Andrew D King, Mayssam Mohammadi Nevisi, Jeremy P Hilton, and Catherine C McGeoch.  Quantum annealing amid local ruggedness and global frustration.  Journal of the Physical Society of Japan, 88(6):061007, 2019.
  * [27] Vladimir Kolmogorov.  Convergent tree-reweighted message passing for energy minimization.  In International Workshop on Artificial Intelligence and Statistics, pages 182–189. PMLR, 2005.
  * [28] Vladimir Kolmogorov and Carsten Rother.  Comparison of energy minimization algorithms for highly connected graphs.  In European Conference on Computer Vision, pages 1–15. Springer, 2006.
  * [29] Nikos Komodakis and Georgios Tziritas.  Approximate labeling via graph cuts based on linear programming.  IEEE transactions on pattern analysis and machine intelligence, 29(8):1436–1453, 2007.
  * [30] Yaroslav Koshka and M. A. Novotny.  2000 Qubit D-Wave Quantum Computer Replacing MCMC for RBM Image Reconstruction and Classification.  In 2018 International Joint Conference on Neural Networks (IJCNN), pages 1–8, July 2018.  ISSN: 2161-4407.
  * [31] Yaroslav Koshka, Dilina Perera, Spencer Hall, and M. A. Novotny.  Determination of the Lowest-Energy States for the Model Distribution of Trained Restricted Boltzmann Machines Using a 1000 Qubit D-Wave 2X Quantum Computer.  Neural Computation, 29(7):1815–1837, 2017.
  * [32] Yaroslav Koshka, Dilina Perera, Spencer Hall, and M.A. Novotny.  Empirical investigation of the low temperature energy function of the Restricted Boltzmann Machine using a 1000 qubit D-Wave 2X.  In 2016 International Joint Conference on Neural Networks (IJCNN), pages 1948–1954, July 2016.  ISSN: 2161-4407.
  * [33] Jiankun Li, Peisen Wang, Pengfei Xiong, Tao Cai, Ziwei Yan, Lei Yang, Jiangyu Liu, Haoqiang Fan, and Shuaicheng Liu.  Practical stereo matching via cascaded recurrent network with adaptive correlation.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 16263–16272, 2022.
  * [34] Junde Li and Swaroop Ghosh.  Quantum-soft QUBO suppression for accurate object detection.  In European Conference on Computer Vision, pages 158–173. Springer, 2020.
  * [35] S. Z. Li.  Markov Random Field Modeling in Computer Vision.  Springer-Verlag, Berlin, Heidelberg, 1995.
  * [36] Andrew Lucas.  Ising formulations of many np problems.  Frontiers in physics, 2:5, 2014.
  * [37] Ayoub Jafari Malekabadi, Mehdi Khojastehpour, and Bagher Emadi.  Disparity map computation of tree using stereo vision system and effects of canopy shapes and foliage density.  Computers and electronics in agriculture, 156:627–644, 2019.
  * [38] Catherine C McGeoch.  Adiabatic quantum computation and quantum annealing: Theory and practice.  Synthesis Lectures on Quantum Computing, 5(2):1–93, 2014.
  * [39] Daniel Scharstein and Richard Szeliski.  A taxonomy and evaluation of dense two-frame stereo correspondence algorithms.  International journal of computer vision, 47(1):7–42, 2002.
  * [40] Jian Sun, Nan-Ning Zheng, and Heung-Yeung Shum.  Stereo matching using belief propagation.  IEEE Transactions on pattern analysis and machine intelligence, 25(7):787–800, 2003.
  * [41] D-Wave Systems.  D-Wave hybrid solver service + advantage: Technology update.  [Online]. Available from: https://www.dwavesys.com/media/m2xbmlhs/14-1048a-a_d-wave_hybrid_solver_service_plus_advantage_technology_update.pdf, June 2022.
  * [42] D-Wave Systems.  Discrete quadratic models.  [Online]. Available from: https://docs.ocean.dwavesys.com/en/latest/concepts/dqm.html, June 2023.
  * [43] Richard Szeliski, Ramin Zabih, Daniel Scharstein, Olga Veksler, Vladimir Kolmogorov, Aseem Agarwala, Marshall Tappen, and Carsten Rother.  A comparative study of energy minimization methods for markov random fields.  In European conference on computer vision, pages 16–29. Springer, 2006.
  * [44] Richard Szeliski, Ramin Zabih, Daniel Scharstein, Olga Veksler, Vladimir Kolmogorov, Aseem Agarwala, Marshall Tappen, and Carsten Rother.  A comparative study of energy minimization methods for markov random fields with smoothness-based priors.  IEEE transactions on pattern analysis and machine intelligence, 30(6):1068–1080, 2008.
  * [45] Marshall F Tappen and William T Freeman.  Comparison of graph cuts with belief propagation for stereo, using identical mrf parameters.  In Computer Vision, IEEE International Conference on, volume 3, pages 900–900. IEEE Computer Society, 2003.
  * [46] Olga Veksler.  Efficient graph-based energy minimization methods in computer vision.  Cornell University, 1999.
  * [47] Middlebury Stereo Vision.  2001 middleburry stereo datasets.  [Online]. Available from: https://vision.middlebury.edu/stereo/data/scenes2001/, June 2022.
  * [48] Middlebury Stereo Vision.  Middleburry MRF implementations.  [Online]. Available from: https://vision.middlebury.edu/MRF/, June 2022.
  * [49] Athanasios Voulodimos, Nikolaos Doulamis, Anastasios Doulamis, Eftychios Protopapadakis, et al.  Deep learning for computer vision: A brief review.  Computational intelligence and neuroscience, 2018, 2018.
  * [50] Martin Wainwright, Tommi Jaakkola, and Alan Willsky.  Tree consistency and bounds on the performance of the max-product algorithm and its generalizations.  Statistics and computing, 14(2):143–166, 2004.
  * [51] Martin J Wainwright, Tommi S Jaakkola, and Alan S Willsky.  Map estimation via agreement on trees: message-passing and linear programming.  IEEE transactions on information theory, 51(11):3697–3717, 2005\.
  * [52] Chaohui Wang, Nikos Komodakis, and Nikos Paragios.  Markov Random Field modeling, inference & learning in computer vision & image understanding: A survey.  Computer Vision and Image Understanding, 117(11):1610–1627, 2013\.
  * [53] Ran Yaacoby, Nathan Schaar, Leon Kellerhals, Oren Raz, Danny Hermelin, and Rami Pugatch.  A comparison between D-Wave and a classical approximation algorithm and a heuristic for computing the ground state of an ising spin glass.  arXiv preprint arXiv:2105.00537, 2021.
  * [54] Tianli Yu, Ruei-Sung Lin, Boaz Super, and Bei Tang.  Efficient message representations for belief propagation.  In 2007 IEEE 11th International Conference on Computer Vision, pages 1–8. IEEE, 2007.
  * [55] Alp Yurtsever, Tolga Birdal, and Vladislav Golyanik.  Q-fw: A hybrid classical-quantum frank-wolfe for quadratic binary optimization.  In European Conference on Computer Vision, pages 352–369. Springer, 2022.
  * [56] Jan-Nico Zaech, Alexander Liniger, Martin Danelljan, Dengxin Dai, and Luc Van Gool.  Adiabatic quantum computing for multi object tracking.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8811–8822, 2022.

##  7 Supplementary Materials: Example*

This part supplements our main paper by providing a simple example to show how
(5) can be modeled and minimized via Quantum Annealing. (4(a)) and (4(a)) show
a pair of $(3\times 4)$-sized stereo images with $D=\\{0,1\\}$. The intensity
values for the left and right images are shown on the pixels. The
corresponding pixel coordinates are illustrated in (4(c)). Without loss of
generality, we ignore the first column of pixels in the left image since
$d_{max}$ is $1$, and we would obtain negative coordinates to match this
column in the right image. The main goal is to compute the disparity map
allocated to the shown red square in (4(a)). Eq. 4(d) shows the ground truth
disparity map.

(a)

(b)

(c)

(d)

Figure 4: (a) the left stereo image, (b) the right stereo image, (c) the pixel
coordinates, (d) the corresponding disparity map.

Considering (1) and (2), we first define $P$ and $N$ follows:

$\displaystyle P=$
$\displaystyle\\{(1,0),(2,0),(3,0),(1,1),(2,1),(3,1),(1,2),(2,2),(3,2)\\}$
$\displaystyle N=$
$\displaystyle\\{\\{(1,0),(2,0)\\},\\{(1,0),(1,1)\\},\\{(2,0),(3,0)\\},$
$\displaystyle\\{(2,0),(2,1)\\},\\{(3,0),(3,1)\\},\\{(1,1),(2,1)\\},$
$\displaystyle\\{(1,1),(1,2)\\},\\{(2,1),(3,1)\\},\\{(2,1),(2,2)\\},$
$\displaystyle\\{(3,1),(3,2)\\},\\{(1,2),(2,2)\\},\\{(2,2),(3,2)\\}\\}.$

The numbers of pixels and disparities are $9$ and $2$, respectively.
Therefore, We define a vector of $18$ binary variables as
$\mathbf{x}=\\{0,1\\}^{18}$:

$\displaystyle\mathbf{x}=($ $\displaystyle
x_{1,0,0},x_{1,0,1},x_{2,0,0},x_{2,0,1},x_{3,0,0},x_{3,0,1},x_{1,1,0},$
$\displaystyle
x_{1,1,1},x_{2,1,0},x_{2,1,1},x_{3,1,0},x_{3,1,1},x_{1,2,0},x_{1,2,1},$
$\displaystyle x_{2,2,0},x_{2,2,1},x_{3,2,0},x_{3,2,1}).$

We set $\lambda=10$ and $\alpha=200$ by which we have
$\alpha>\left(\sum_{(i,j)\in P}\max\\{\theta_{i,j}(d)|d\in
D\\}\right)+\lambda|N|$. The QUBO model (5) is formulated as follows:

$\displaystyle H(\mathbf{x})$ $\displaystyle=200\sum_{(i,j)\in
P}\left(1-\sum_{d\in D}x_{i,j,d}\right)^{2}+\sum_{(i,j)\in P}\sum_{d\in
D}\theta_{i,j}(d)x_{i,j,d}$
$\displaystyle+10\sum_{\\{(i,j),(i^{\prime},j^{\prime})\\}\in N}\sum_{d_{1}\in
D}\sum_{d_{2}\in
D}\delta(d_{1},d_{2})x_{i,j,d_{1}}x_{i^{\prime},j^{\prime},d_{2}}.$

The QUBO model $H(\mathbf{x})$ has three terms denoted by $H_{1}$, $H_{2}$,
and $H_{3}$ from left to right, respectively. The following shows each term’s
expansion separately. We then add them all at the end. We start with the first
term denoted by $H_{1}$.

$\displaystyle H_{1}(\mathbf{x})=200($ $\displaystyle-
x_{1,0,0}-x_{1,0,1}+2x_{1,0,0}x_{1,0,1}+1$ $\displaystyle-
x_{2,0,0}-x_{2,0,1}+2x_{2,0,0}x_{2,0,1}+1$ $\displaystyle-
x_{3,0,0}-x_{3,0,1}+2x_{3,0,0}x_{3,0,1}+1$ $\displaystyle-
x_{1,1,0}-x_{1,1,1}+2x_{1,1,0}x_{1,1,1}+1$ $\displaystyle-
x_{2,1,0}-x_{2,1,1}+2x_{2,1,0}x_{2,1,1}+1$ $\displaystyle-
x_{3,1,0}-x_{3,1,1}+2x_{3,1,0}x_{3,1,1}+1$ $\displaystyle-
x_{1,2,0}-x_{1,2,1}+2x_{1,2,0}x_{1,2,1}+1$ $\displaystyle-
x_{2,2,0}-x_{2,2,1}+2x_{2,2,0}x_{2,2,1}+1$ $\displaystyle-
x_{3,2,0}-x_{3,2,1}+2x_{3,2,0}x_{3,2,1}+1)$

Next, we expand the second term as $H_{2}$.

$\displaystyle H_{2}(\mathbf{x})$
$\displaystyle=|I_{l}(1,0)-I_{r}(1,0)|x_{1,0,0}+|I_{l}(1,0)-I_{r}(0,0)|x_{1,0,1}$
$\displaystyle+|I_{l}(2,0)-I_{r}(2,0)|x_{2,0,0}+|I_{l}(2,0)-I_{r}(1,0)|x_{2,0,1}$
$\displaystyle+|I_{l}(3,0)-I_{r}(3,0)|x_{3,0,0}+|I_{l}(3,0)-I_{r}(2,0)|x_{3,0,1}$
$\displaystyle+|I_{l}(1,1)-I_{r}(1,1)|x_{1,1,0}+|I_{l}(1,1)-I_{r}(0,1)|x_{1,1,1}$
$\displaystyle+|I_{l}(2,1)-I_{r}(2,1)|x_{2,1,0}+|I_{l}(2,1)-I_{r}(1,1)|x_{2,1,1}$
$\displaystyle+|I_{l}(3,1)-I_{r}(3,1)|x_{3,1,0}+|I_{l}(3,1)-I_{r}(2,1)|x_{3,1,1}$
$\displaystyle+|I_{l}(1,2)-I_{r}(1,2)|x_{1,2,0}+|I_{l}(1,2)-I_{r}(0,2)|x_{1,2,1}$
$\displaystyle+|I_{l}(2,2)-I_{r}(2,2)|x_{2,2,0}+|I_{l}(2,2)-I_{r}(1,2)|x_{2,2,1}$
$\displaystyle+|I_{l}(3,2)-I_{r}(3,2)|x_{3,2,0}+|I_{l}(3,2)-I_{r}(2,2)|x_{3,2,1}$
$\displaystyle H_{2}(\mathbf{x})$
$\displaystyle=50x_{1,0,0}+50x_{2,0,1}+50x_{2,1,0}+50x_{3,1,1}+50x_{1,2,0}+50x_{2,2,1}.$

Finally, we compute the third term as $H_{3}$:

$\displaystyle H_{3}(\mathbf{x})$
$\displaystyle=10(x_{1,0,0}x_{2,0,1}+x_{1,0,1}x_{2,0,0}$
$\displaystyle+x_{1,0,0}x_{1,1,1}+x_{1,0,1}x_{1,1,0}+x_{2,0,0}x_{3,0,1}$
$\displaystyle+x_{2,0,1}x_{3,0,0}+x_{2,0,0}x_{2,1,1}+x_{2,0,1}x_{2,1,0}$
$\displaystyle+x_{3,0,0}x_{3,1,1}+x_{3,0,1}x_{3,1,0}+x_{1,1,0}x_{2,1,1}$
$\displaystyle+x_{1,1,1}x_{2,1,0}+x_{1,1,0}x_{1,2,1}+x_{1,1,1}x_{1,2,0}$
$\displaystyle+x_{2,1,0}x_{3,1,1}+x_{2,1,1}x_{3,1,0}+x_{2,1,0}x_{2,2,1}$
$\displaystyle+x_{2,1,1}x_{2,2,0}+x_{3,1,0}x_{3,2,1}+x_{3,1,1}x_{3,2,0}$
$\displaystyle+x_{1,2,0}x_{2,2,1}+x_{1,2,1}x_{2,2,0}+x_{2,2,0}x_{3,2,1}$
$\displaystyle+x_{2,2,1}x_{3,2,0}).$

Adding the three terms together, we have the main QUBO model as follows:

$\displaystyle H(\mathbf{x})=$
$\displaystyle-150x_{1,0,0}-200x_{1,0,1}-200x_{2,0,0}-150x_{2,0,1}$
$\displaystyle-200x_{3,0,0}-200x_{3,0,1}-200x_{1,1,0}-200x_{1,1,1}$
$\displaystyle-150x_{2,1,0}-200x_{2,1,1}-200x_{3,1,0}-150x_{3,1,1}$
$\displaystyle-150x_{1,2,0}-200x_{1,2,1}-200x_{2,2,0}-150x_{2,2,1}$
$\displaystyle-200x_{3,2,0}-200x_{3,2,1}$
$\displaystyle+400x_{1,0,0}x_{1,0,1}+400x_{2,0,0}x_{2,0,1}$
$\displaystyle+400x_{3,0,0}x_{3,0,1}+400x_{1,1,0}x_{1,1,1}$
$\displaystyle+400x_{2,1,0}x_{2,1,1}+400x_{3,1,0}x_{3,1,1}$
$\displaystyle+400x_{1,2,0}x_{1,2,1}+400x_{2,2,0}x_{2,2,1}$
$\displaystyle+400x_{3,2,0}x_{3,2,1}+1800$
$\displaystyle+10x_{1,0,0}x_{2,0,1}+10x_{1,0,1}x_{2,0,0}+10x_{1,0,0}x_{1,1,1}$
$\displaystyle+10x_{1,0,1}x_{1,1,0}+10x_{2,0,0}x_{3,0,1}+10x_{2,0,1}x_{3,0,0}$
$\displaystyle+10x_{2,0,0}x_{2,1,1}+10x_{2,0,1}x_{2,1,0}+10x_{3,0,0}x_{3,1,1}$
$\displaystyle+10x_{3,0,1}x_{3,1,0}+10x_{1,1,0}x_{2,1,1}+10x_{1,1,1}x_{2,1,0}$
$\displaystyle+10x_{1,1,0}x_{1,2,1}+10x_{1,1,1}x_{1,2,0}+10x_{2,1,0}x_{3,1,1}$
$\displaystyle+10x_{2,1,1}x_{3,1,0}+10x_{2,1,0}x_{2,2,1}+10x_{2,1,1}x_{2,2,0}$
$\displaystyle+10x_{3,1,0}x_{3,2,1}+10x_{3,1,1}x_{3,2,0}+10x_{1,2,0}x_{2,2,1}$
$\displaystyle+10x_{1,2,1}x_{2,2,0}+10x_{2,2,0}x_{3,2,1}+10x_{2,2,1}x_{3,2,0}.$

Giving $H(\mathbf{x})$ to the D-Wave Ocean SDK for the QPU minimization, we
obtain the optimal solution
$\mathbf{x}^{*}=\operatorname*{arg\,min}_{\mathbf{x}}H(\mathbf{x})$ as
follows:

  * •

$x^{*}_{1,0,0}=0,$

  * •

$x^{*}_{1,0,1}=1,$

  * •

$x^{*}_{2,0,0}=1,$

  * •

$x^{*}_{2,0,1}=0,$

  * •

$x^{*}_{3,0,0}=1,$

  * •

$x^{*}_{3,0,1}=0,$

  * •

$x^{*}_{1,1,0}=0,$

  * •

$x^{*}_{1,1,1}=1,$

  * •

$x^{*}_{2,1,0}=0,$

  * •

$x^{*}_{2,1,1}=1,$

  * •

$x^{*}_{3,1,0}=1,$

  * •

$x^{*}_{3,1,1}=0,$

  * •

$x^{*}_{1,2,0}=0,$

  * •

$x^{*}_{1,2,1}=1,$

  * •

$x^{*}_{2,2,0}=1,$

  * •

$x^{*}_{2,2,1}=0,$

  * •

$x^{*}_{3,2,0}=1,$

  * •

$x^{*}_{3,2,1}=0,$

which match the corresponding ground-truth disparities shown in (4(d)). We
used the D-Wave default parameter settings for the hardware properties and
initialized the number of sample-reads as $1000$. Given a pixel $(i,j)\in P$,
if $x^{*}_{i,j,d}=1$ for $d\in D$, then $d$ is the allocated disparity to the
pixel $(i,j)$. Therefore, we have the following disparities for the pixels:

  * •

$(1,0)\leftarrow 1$

  * •

$(2,0)\leftarrow 0$

  * •

$(3,0)\leftarrow 0$

  * •

$(1,1)\leftarrow 1$

  * •

$(2,1)\leftarrow 1$

  * •

$(3,1)\leftarrow 0$

  * •

$(1,2)\leftarrow 1$

  * •

$(2,2)\leftarrow 0$

  * •

$(3,2)\leftarrow 0$

Eq. 5 illustrates the corresponding D-Wave minor embedding for the defined
QUBO, obtained by the D-Wave Inspector tool. Furthermore, we provide the
timing information in (3).

Table 3: The D-Wave timing information for the given QUBO example Type | Time ($\mu s$)
---|---
$qpu\\_sampling\\_time$ | $104760.00$
$qpu\\_anneal\\_time\\_per\\_sample$ | $20.00$
$qpu\\_readout\\_time\\_per\\_sample$ | $64.22$
$qpu\\_access\\_time$ | $120523.57$
$qpu\\_access\\_overhead\\_time$ | $4687.43$
$qpu\\_programming\\_time$ | $15763.57$
$qpu\\_delay\\_time\\_per\\_sample$ | $20.54$
$total\\_post\\_processing\\_time$ | $852.00$
$post\\_processing\\_overhead\\_time$ | $852.00$

(a)

(b)

Figure 5: D-Wave minor embedding for the given QUBO example. (a) the QUBO
graph, and (b) the QPU graph.