Best Subset Selection: Statistical Computing Meets Quantum Computing

To facilitate the discussion in the paper, we review some essential notations and definitions used in quantum computing that may be alien to the audience of statistics. We refer to Nielsen and Chuang, (2010) and Schuld and Petruccione, (2018) for more detailed tutorials of quantum computing. We start by introducing the linear algebra notations used in quantum mechanics, which were invented by Paul Dirac. In Dirac’s notation, a column vector ${\bf a}$ is denoted by $\ket{a}$ which reads as ${\bf a}$ in a ‘ket’. Typical vector space of interest in quantum computing is a Hilbert space $\mathcal{H}$ of dimension $D=2^{p}$ with a positive integer $p$. For a vector $\ket{a}\in\mathcal{H}$, we denote its dual vector as $\bra{a}$ (reads as $\bf a$ in a ‘bra’), which is an element of $\mathcal{H}^{*}$. Besides, $\mathcal{H}$ and $\mathcal{H}^{*}$ together naturally induce an inner product $\braket{a}{b}=\braket{{\bf a},{\bf b}}$, which is also known as a ‘bra-ket’. We say $\ket{a}$ is a unit vector if $\braket{a}{a}=1$.

The framework of quantum computing resides in a state-space postulate which describes the state of a system by a unit vector in a Hilbert space. A state of a qubit can be described by a unit vector $\ket{\psi}$ in a two-dimensional Hilbert space. A superposition, as the linear combination of the states in one qubit, can be defined as

where coefficients $\phi_{0}$ and $\phi_{1}$ satisfy $|\phi_{0}|^{2}+|\phi_{1}|^{2}=1$. Now suppose that we have two qubits and the system has four orthonormal states $\ket{00},\ket{01},\ket{10},\ket{11}$. Then, a superposition can be represented as

where coefficients $\phi_{0}$, $\phi_{1}$, $\phi_{2}$, $\phi_{3}$ satisfiy $|\phi_{0}|^{2}+|\phi_{1}|^{2}+|\phi_{2}|^{2}+|\phi_{3}|^{2}=1$. Analogously, a quantum computer of $p$ qubits can represent a state of a system by a unit vector $\ket{\psi}$ in a $D=2^{p}$ dimensional Hilbert space $\mathcal{H}$. Let ${\mathcal{B}}=\{\ket{b_{i}}\}_{i=0}^{D-1}$ be an orthonormal basis of $\mathcal{H}$. Every superposition $\ket{\psi}\in\mathcal{H}$ can be represented as

where $\phi_{0},\ \ldots,\ \phi_{D-1}$ is a set of coefficients with $\phi_{i}=\braket{b_{i}}{\psi}$ and $\sum_{i=0}^{D-1}|\phi_{i}|^{2}=1$.

Another salient feature of quantum computing is the measurement of a quantum system yields a probabilistic outcome rather than a deterministic one. Suppose the superposition $\ket{\psi}$ in (1) is measured, it collapses to a random state in ${\mathcal{B}}$. In addition, the probability we observe $\ket{b_{i}}$ is $|\phi_{i}|^{2}$ for $i=0,\ldots,D-1$.

Different from a classical computer with logic gates, e.g., AND and NOT gates, a quantum computer is built with quantum gates. These quantum gates not only have the logic gates mentioned above but also include the new gates that are not available in classical computers. For example, the Hadamard gate enables the rotations and reflections of the input states. Some quantum algorithms that make use of these new gates have shown some remarkable results that significantly outperform the algorithms in classical computers.

The best subset selection problem can be decomposed into two subproblems: (1) What is the criterion for the best? (2) Given the criterion, which subset is the best? For the first subproblem, there are many criteria, e.g, AIC, BIC, cross-validation, and mean squared error. For the second subproblem, searching for the best is challenging if the number of candidate subsets is super-large. We shall review Grover’s algorithm, which is the most popular quantum searching algorithm for the second subproblem.

Recall that given $p$ predictors, there are $D=2^{p}$ candidate subsets, from which we aim to find the best subset. Once the best criterion is chosen and applied to each subset, we have $D=2^{p}$ criterion values computed for all subsets. We aim to find the solution (e.g. the smallest or the largest value) from $D=2^{p}$ elements, which are indexed by $0,1,\ldots,D-1$. To ease the presentation, we assume there are no tied solutions although the discussions in this paper can be easily generalized to such cases. In the language of quantum computing, each index is regarded as a state, which can be modeled by a unit vector in an orthonormal basis ${\cal D}=\{\ket{i}\})_{i=0}^{D-1}$ of a $D$ dimensional Hilbert space $\mathcal{H}$.

Best subset selection has been a classical and fundamental method for linear regressions (e.g. Hocking and Leslie,, 1967; Beale et al.,, 1967; Shen et al.,, 2012; Fan et al.,, 2020). When the covariates contain redundant information for the response, selecting a parsimonious best subset of covariates may improve the model interpretability as well as the prediction performance. The best subset selection for linear regression can be formulated as

where $\bm{Y}=(y_{1},\ \ldots,\ y_{n})^{\mathrm{\scriptstyle T}}\in\mathbb{R}^{n}$ is a response vector, $\bm{X}=(x_{i,j})_{1\leq i\leq n,1\leq j\leq p}\in\mathbb{R}^{n\times p}$ is a design matrix of covariates and $\bm{\beta}=(\beta_{1},\ \ldots,\ \beta_{p})^{\mathrm{\scriptstyle T}}\in\mathbb{R}^{p}$ is a vector of unknown regression coefficients. Throughout this paper, we assume $\bm{Y}$ and $\bm{X}$ are centered for the simplicity of presentation. Further, $\|\bm{\beta}\|_{0}$ is the $\ell_{0}$ norm of $\bm{\beta}$ and $0\leq t\leq\min{(n,p)}$ is a subset size. In general, solving (4) is a non-convex optimization problem. Without a convex relaxation or a stage-wise approximation, the optimization involves a combinatorial search over all subsets of sizes $\{1,\ \ldots,\ \min(n,p)\}$ and hence is an NP-hard problem (Welch,, 1982; Natarajan,, 1995). Therefore, seeking the exact solution of (4) is usually computationally intractable when $p$ is moderate or large.

To overcome the computational bottleneck of the best subset selection, many alternative and approximate approaches have been developed and carefully studied. Stepwise regression methods (e.g. Efroymson,, 1960; Draper and Smith,, 1966) sequentially add or eliminate covariates according to their marginal contributions to the response condition on the existing model. Stepwise regression methods speed up the best subset selection as they only optimize (4) over a subset of all candidates, and hence their solutions may not coincide with the solution of (4). However, stepwise regression methods suffer from the instability issue in high dimensional regime (Breiman,, 1996). Regularized regression methods suggest to relax the nonconvex $\ell_{0}$ norm in (4) by various convex or nonconcave alternatives. A celebrated member in the house is LASSO (Tibshirani,, 1996) which solves an $\ell_{1}$ regularized problem instead. Elegant statistical properties and promising numerical performance soon made regularized regression a popular research topic in statistics, machine learning and other data science related areas. We refer to basis pursuit (Chen and Donoho,, 1994), SCAD (Fan and Li,, 2001), elastic net (Zou and Hastie,, 2005), adaptive LASSO (Zou,, 2006), Dantzig selector (Candes and Tao,, 2007), nonnegative garrotte (Yuan and Lin,, 2007), and MCP (Zhang,, 2010), among many others. Though regularized regression methods can perform equivalent to or even better than the solution of (4) in certain scenarios (e.g. Zhao and Yu,, 2006; Fan et al.,, 2014; Hastie et al.,, 2020), they are not universal replacements for the best subset selection method. Nevertheless, solving (4) efficiently remains a statistically attractive and computationally challenging problem.

Recently, Bertsimas et al., (2016) has reviewed the best subset selection problem from a modern optimization point of view. Utilizing the algorithmic advances in mixed integer optimization, Bertsimas et al., (2016) proposed a highly optimized solver of (4) which is scalable to relatively large $n$ and $p$. Later, Hazimeh and Mazumder, (2018) considered a new hierarchy of necessary optimality conditions and propose to solve (4) by adding extra convex regularizations. Then, Hazimeh and Mazumder, (2018) developed an efficient algorithm based on coordinate descent and local combinatorial optimization. Such recent optimization developments push the frontier of computation for the best subset selection problems by a big margin.

Let us consider the best subset selection problem (4) and assume $n\geq p$ such that all $D=2^{p}$ subsets of $p$ covariates are attainable. When $n<p$, we only need to search over $\sum_{t=0}^{n}{p\choose t}<2^{p}$ subsets which is a less challenging combinatorial search problem.

Intuitively, one can think of QAS as a “quantum elevator” that starts at a random floor of a high tower and aims to descent to the ground floor. In each operation, a quantum machine will randomly decide if this elevator stays at the current floor or goes down to a random lower-level floor. Besides, the probability of going down will gradually increase as the number of operations increases. Thus, it is not hard to imagine that the “quantum elevator” can descent to the ground floor with a high success probability after a large enough amount of operations. Our intuition of QAS is justified by the following theorem. Theorem 3.1 states that within $O(\log_{2}D)$ iterations, QAS finds the sole solution state for any pre-specified success probability greater than 50%.

Suppose that, in the $m$th iteration of Algorithm 2, the current benchmark state is $\ket{w_{m}}$ with $g(\ket{w_{m}}$ being the $r_{m}$th smallest among $\{g(\ket{i})\}_{i=0}^{D-1}$. Theorem 3.2 below shows the expected number of Grover’s operations that Algorithm 2 needs to update $g(\ket{w_{m}})$ is on the order $O(\sqrt{D/r_{m}})$. This together with Theorem 3.1 imply the computational complexity upper bound of QAS is of order $O(\sqrt{D}\log_{2}D)$ which is only a $\log_{2}D$ factor larger than the theoretical lower bound of any oracular quantum search algorithm (Bennett et al.,, 1997). Notice that, QAS achieves a near optimal computational efficiency as oracular quantum search algorithms without using any oracle information.

The success of QAS relies on the efficient comparison of two states through the state loss function defined in (5), which is the MSE of linear regression. Recently, a line of studies (Wiebe et al.,, 2012; Rebentrost et al.,, 2014; Wang,, 2017; Zhao et al.,, 2019) focuses on formulating linear regression as a matrix inversion problem, which can be tackled by the quantum algorithm to solve linear systems of equations (Harrow et al.,, 2009). Their goal is to find a series of quantum states whose amplitudes represent the ordinary least squares estimator of linear regression coefficients. Although their algorithms can encode such quantum states efficiently, “it may be exponentially expensive to learn via tomography” (Harrow et al.,, 2009). As a result, most existing quantum linear regression algorithms require the input design matrix to be sparse. Otherwise, the computation is slow and the estimation accuracy may suffer from noise accumulation.

In this paper, we investigate a novel quantum linear prediction approach that is based on the singular-value decomposition. The proposed method avoids the matrix inversion problem. Instead, we estimate the inverse singular values by utilizing a recently developed quantum tomography technique (Lloyd et al.,, 2014). Let $\{j_{1},\ \ldots,\ j_{d}\}$ be a subset of $\{1,\ \ldots,\ p\}$ with $d\leq n$. Denote $\bm{X}_{d}\in\mathbb{R}^{n\times d}$ the sub-matrix of $\bm{X}$ corresponding to the subset $\{j_{1},\ \ldots,\ j_{d}\}$. Let $\widetilde{\bm{x}}$ be a new observation of the $d$ selected covariates. Then, the corresponding response value $\widetilde{y}$ can be predicted by

Suppose that the rank of $\bm{X}_{d}$ is $R\leq d$, a compact singular value decomposition of $\bm{X}_{d}$ yields

where $\bm{\Sigma}=\mathrm{diag}\{\sigma_{1},\ \ldots,\ \sigma_{R}\}$ is a diagonal matrix of $R$ positive singular values, and $\bm{U}=(\bm{u}_{1},\ \ldots,\ \bm{u}_{R})\in\mathbb{R}^{n\times R}$ and $\bm{V}=(\bm{v}_{1},\ \ldots,\ \bm{v}_{R})\in\mathbb{R}^{d\times R}$ are matrices of left and right singular vectors, respectively. Further, we have $\bm{U}\bm{U}^{\mathrm{\scriptstyle T}}=\bm{V}\bm{V}^{\mathrm{\scriptstyle T}}=\bm{I}_{R}$. Then, the linear predictor in (7) can be represented by

Motivated by the quantum principle component analysis (QPCA) (Lloyd et al.,, 2014) and the inverted singular value problem (Schuld et al.,, 2016), we propose to calculate (8) by a quantum linear prediction (QLP) procedure, which can be summarized as the following three steps.

We encode $\bm{X}_{d}$, $\bm{v}_{r}$ and $\bm{u}_{r}$, $r=1,\ \ldots,\ R$, into a quantum system as

where $\sum_{i}\sum_{j}|a_{i,j}|^{2}=\sum_{i}|\mu_{r,i}|^{2}=\sum_{j}|\nu_{r,j}|^{2}=1$, and $\{\ket{i}\}_{i=1}^{n}$ and $\{\ket{j}\}_{j=1}^{d}$ are orthonormal quantum bases representing the linear space spanned by the left and right singular vectors of $\bm{X}_{d}$, respectively.

Then, we encode $\bm{X}_{d}^{\mathrm{\scriptstyle T}}\bm{X}_{d}$ and the compact singular value decomposition of $\bm{X}_{d}$ into a quantum system as

Similarly, we can encode $\bm{y}$ and $\widetilde{\bm{x}}$ into a quantum system as

with $\sum_{\eta}|b_{\eta}|^{2}=\sum_{\gamma}|c_{\gamma}|^{2}=1$, and $\{\ket{\eta}\}_{\eta=1}^{n}$ and $\{\ket{\gamma}\}_{\gamma=1}^{d}$ are orthonormal quantum bases.

To calculate (8) with a quantum computer, we first need to calculate the inverse of singular values of $\bm{X}_{d}$. According to the ideas of QPCA (Lloyd et al.,, 2014), we apply ${\mathbf{\rho}}$ to $\ket{{\bm{X}_{d}}}$ and follow the quantum phase estimation algorithm (e.g. Shor,, 1994; Szegedy,, 2004; Wocjan et al.,, 2009) which yields

where $\lambda_{r}=\sigma_{r}^{2}$, $r=1,\ \ldots,\ R$, are the eigenvalues of ${\mathbf{\rho}}$ which is encoded in a quantum register $\ket{\lambda_{r}}$. Then, by adding an extra qubit and rotating (12) condition on the eigenvalue register $\ket{\lambda_{r}}$, we have

where $c$ is a constant to ensure the inverse of eigenvalues are no larger than 1.

We can repeatedly perform a conditional measurement of (13) until it is in state $\ket{1}$. After that, we can uncompute and discard the eigenvalue register which results in

where $p(1)=\sum_{r=1}^{R}|\frac{c}{\lambda_{r}}|^{2}$ is the probability of the qubit in (13) collapsing to state $\ket{1}$.

Follow the notations in (9) –(11), we can rewrite (8) as

Motivated by the strategy in Schuld et al., (2016), we write (15) into some entries of a qubit so that it can be assessed by a single measurement. To this end, we construct two states $\ket{\psi_{1}}$ and $\ket{\psi_{2}}$ that are entangled with a qubit, i.e.

where $\ket{\psi_{1}}$ is defined in (14) and $\ket{\psi_{2}}=\ket{{\bm{y}}}\ket{{\widetilde{\bm{x}}}}$. Then, the off-diagonal elements of this qubit’s density matrix read

which contain the desired results (8) up to a known normalization factor.

Given a testing set of $n_{t}$ observations, i.e. $\{\widetilde{y}_{i},\widetilde{\bm{x}}_{i}\}_{i=1}^{n_{t}}$, the prediction error can be calculated as

where $\ket{i}$ is a quantum state represents $\{x_{1},\ \ldots\ ,x_{d}\}$ and $\widehat{y}_{i}$ is the predictor of $\widetilde{y}_{i}$ which can be calculated follow the QLP procedure introduced above. In addition, the computational complexity of running QLP over a training set of size $n$ and a testing set of $n_{t}$ is approximately of the order $O\left((n+n_{t})\log_{2}d\right)$ which is faster than the classical computing algorithm that is usually of order $O\left((n+n_{t})d\right)$.

The proposed hybrid quantum-classical strategy involves two tuning parameters: the number of quantum nodes $K$; and a learning rate $\lambda\in(0,1)$ which will be passed to each quantum node to implement Algorithm 2. As suggested by Theorem 4.1, the success probability of Algorithm 3 converges to 1 as $K$ increases when the success probability of each quantum node surpasses 0.5. Therefore, one should choose $K$ as large as possible subject to the availability of quantum computing resources. On the other hand, the learning rate $\lambda$ controls the “step size" of each iteration in Algorithm 2. In this subsection, we use an experiment to assess the performance and the sensitivity of Algorithm 3 with respect to the choices of $K$ and $\lambda$. To focus on the selection of tuning parameters, we skip the stage 1 in Algorithm 3. Instead, we generate $\bm{G}_{D}$ as an i.i.d. sample of size $D=32$ from the uniform distribution $U[0,1]$. Then we use the stages 2 and 3 in Algorithm 3 to find the minimum element in $\bm{G}_{D}$. We set $K=1$, $3$ or $5$, and let $\lambda$ be a sequence of grid points between 0.40 and 0.60 with a step size of 0.01. For each pair of $K$ and $\lambda$, we simulate 200 replications. The performance is measured by the accuracy rate which is defined as

The results are presented in Figure 5. The accuracy rates for $K=1,3$ or $5$ are plotted as blue, orange, and green dots with different symbols. The solid lines are corresponding smoothed curves. According to Figure 5, we find that increasing the number of quantum nodes $K$ can improve the accuracy rate though the improvement is marginal when $K$ is large enough, e.g. $K=5$. Besides, we observe that the highest accuracy rates are attained with some $\lambda$ between 0.5 and 0.55 for both choices of $K$. One can notice that the accuracy rates are close to 1 when $\lambda\in(0.5,0.55)$ and $K\geq 3$ which can be considered as a rule-of-thumb recommendation. Further, the smoothed curves in Figure 5 indicate that Algorithm 3 is not very sensitive for the selection of tuning parameters.

In this subsection, we assess the best subset selection performance of the proposed hybrid quantum-classical strategy. We consider a linear regression model

where $y_{i}\in\mathbb{R}$, $\bm{x}_{i}\in\mathbb{R}^{p}$, $\bm{\beta}^{\ast}=(\beta_{1},\ \ldots\ ,\beta_{p})^{\mathrm{\scriptstyle T}}\in\mathbb{R}^{p}$, and $\epsilon_{i}\in\mathbb{R}$. First, we draw $\{\bm{x}_{i}\}_{i=1}^{n}$ as an i.i.d. sample from a multivariate normal distribution $N_{p}({\bf 0},\bm{\Sigma})$, where the $(i,j)$th entry of $\bm{\Sigma}$ equals $\rho^{|i-j|}$ for some $\rho\in(0,1)$. Then, we draw $\{\epsilon_{i}\}_{i=1}^{n}$ as an i.i.d. sample from a normal distribution $N(0,\sigma^{2})$. Notice that, $\rho$ controls the autocorrelation level among covariates, and $\bm{\beta}^{\ast}$, $\rho$ and $\sigma$ together control the signal to noise ratio (SNR), i.e. $\text{SNR}=\bm{\beta}^{\ast\mathrm{\scriptstyle T}}\bm{\Sigma}\bm{\beta}^{\ast}/\sigma^{2}$. Let $s<p$ be a positive integer indicating the model sparsity. We set the first $s$ elements of $\bm{\beta}^{\ast}$ to be 1 and the rest $p-s$ elements to be 0.

Subject to the capacity of the quantum computing system, we set $n=100$, $p=7$, and $s=4$. We choose the autocorrelation level $\rho\in\{0.25,0.5\}$ and the signal to noise ratio $\text{SNR}\in\{0.5,1.0,2.0,3.0\}$, respectively. For each scenario, we simulate 200 replications. We implement the hybrid quantum-classical strategy as proposed in Algorithm 3 with $K=3$ and $\lambda=0.55$ to select the best subset of the linear regression model. We denote this method as QAS since it implements the quantum adaptive search in stage 2. Besides, we consider two variants of Algorithm 3 as two competitors. The first competitor, denoted as Grover Oracle, replaces QAS in stage 2 by Grover’s algorithm with a known oracle state. The second competitor, denoted as Grover Random, replaces QAS in stage 2 by Grover’s algorithm with a randomly selected oracle state. Notice that Grover Oracle is an oracular method and not applicable in practice since the oracle state is not observable.

For each replication, we record the false positives (FP) and false negatives (FN) of each method. FP is the number of inactive covariates that are selected into the model and FN is the number of active covariates that are ignored in the selected model. The box-plots of FP and FN over 200 replications are reported in Figure 6 below. According to Figure 6, Grover Oracle performs the best among the three competitors which is not surprising as it utilizes the oracle information of the true best subset. The proposed QAS method, as a non-oracular approach, performs almost identical to Grover Oracle in nearly all scenarios. QAS is slightly outperformed by Grover Oracle when the correlation level is high ($\rho=0.5$) and the signal to noise ratio is small ($\text{SNR}=0.5$), which is the most challenging scenario. In contrast, Grover Random performs as bad as random guesses in all scenarios. This observation, which is inline with the discussions in Section 2.2.3, shows oracular quantum search algorithms are not suitable to statistical learning problems.