A Quantum-inspired Algorithm for General Minimum Conical Hull Problems

TL;DR

This paper introduces a quantum-inspired classical algorithm that significantly speeds up solving general minimum conical hull problems, especially in high-dimensional data scenarios, by achieving exponential speedup over existing methods.

Contribution

The paper presents a novel sublinear classical algorithm for minimum conical hull problems with exponential speedup over traditional polynomial-time algorithms.

Findings

Achieves exponential speedup over DCA

Runtime is polynomial in rank and polylogarithmic in size

Effective for high-dimensional machine learning tasks

Abstract

A wide range of fundamental machine learning tasks that are addressed by the maximum a posteriori estimation can be reduced to a general minimum conical hull problem. The best-known solution to tackle general minimum conical hull problems is the divide-and-conquer anchoring learning scheme (DCA), whose runtime complexity is polynomial in size. However, big data is pushing these polynomial algorithms to their performance limits. In this paper, we propose a sublinear classical algorithm to tackle general minimum conical hull problems when the input has stored in a sample-based low-overhead data structure. The algorithm's runtime complexity is polynomial in the rank and polylogarithmic in size. The proposed algorithm achieves the exponential speedup over DCA and, therefore, provides advantages for high dimensional problems.

Equations

$$\mathrm{MCH}(\mathbf{Y},\mathbf{X}):=\arg \min \{|\mathcal{A}|:\text{cone}({\mathbf{Y}}_{\mathcal{A}})\supset \text{cone}(\mathbf{X})\},$$

$$\text{MCH}({\mathbf{Y}}_t,{\mathbf{X}}_t):=\arg \min \{|{\mathcal{A}}_t|:cone({\mathbf{Y}}_{{\mathcal{A}}_t})\supset cone({\mathbf{X}}_t)\},$$

$${\hat{g}}_i:=\frac{1}{p}\sum _{t=1}^p\text{mathbbm}{1}_{{\mathcal{A}}_t}(\mathbf{Y}(i;:)),$$

$$\text{MCH}({\mathbf{Y}}_t,{\mathbf{X}}_t):={\mathcal{A}}_t=\left\{\arg \underset{i\in [n_{\text{scalebox}0.49Y}]}{\min }({\mathbf{Y}}_t(i,:)-\underset{j\in [n_{\text{scalebox}0.49X}]}{\max }{\mathbf{X}}_t(j,:){)}_{+}\right\},$$

$$\text{MCH}({\mathbf{Y}}_t,{\mathbf{X}}_t):={\mathcal{A}}_t=\left\{\arg \underset{i\in [n_{\text{scalebox}0.49Y}]}{\min }{\left({\mathcal{P}}_{{\mathbf{Y}}_t}(i)-{\xi }_t\underset{j\in [n_{\text{scalebox}0.49X}]}{\max }({\mathcal{P}}_{{\mathbf{X}}_t}(j))\right)}_{+}\right\},\forall t\in [p],$$

$$\tilde{\mathbf{V}}(:,i)\equiv {\tilde{\mathbf{v}}}^{(i)}:=\frac{\mathbf{R}{\mathbf{\omega }}^{(i)}}{{\sigma }^{(i)}},$$

$$\max \left\{\tilde{\mathcal{O}}(\frac{85^6{k}_{\text{scalebox}0.49X}^6{\kappa }_{\text{scalebox}0.49X}^{12}\Vert \mathbf{X}{\Vert }_F^6}{9^3{ϵ}^6}),\tilde{\mathcal{O}}(\frac{85^6{k}_{\text{scalebox}0.49Y,t}^6{\kappa }_{\text{scalebox}0.49Y}^{12}\Vert \mathbf{Y}{\Vert }_F^6}{9^3{ϵ}^6})\right\}.$$

$$s=\frac{85^2{k}^2{\kappa }^4\ln (8n/\eta )\Vert \mathbf{H}{\Vert }_F^2}{9{ϵ}^2},$$

$$Z(j,l)=\frac{\Vert {\tilde{\mathbf{v}}}^{(i)}\Vert \Vert {\mathbf{B}}_t\Vert \mathbf{H}(j,l)}{{\tilde{\mathbf{v}}}^{(i)}(j){\mathbf{B}}_t(l)}$$

$${\mathcal{P}}_Z(Z(j,l))=\frac{|{\tilde{\mathbf{v}}}^{(i)}(j){|}^2|{\mathbf{B}}_t(l){|}^2}{\Vert {\tilde{\mathbf{v}}}^{(i)}{\Vert }^2\Vert {\mathbf{B}}_t{\Vert }^2}.$$

$$N_Z=N_p\times N_q\sim \mathcal{O}\left(\frac{\Vert \mathbf{H}{\Vert }_F\Vert {\mathbf{B}}_t\Vert \Vert {\tilde{\mathbf{v}}}^{(i)}\Vert }{{ϵ}^2}\log (\frac{1}{\delta })\right).$$

$$|{\bar{z}}^{\ast }-{\tilde{\mathbf{q}}}_t(i)|\le ϵ.$$

$$\mathcal{O}\left(\frac{\Vert \mathbf{H}{\Vert }_F\Vert {\mathbf{B}}_t\Vert \Vert {\tilde{\mathbf{v}}}^{(i)}\Vert }{{ϵ}^2}(L(\mathbf{H})+Q({\mathbf{B}}_t)+Q(\tilde{\mathbf{V}}))\log (\frac{1}{\delta })\right).$$

$$\mathcal{O}\left(\frac{(L(\mathbf{H})+Q({\mathbf{B}}_t)+Q(\tilde{\mathbf{V}}))}{{ϵ}^2}\log (\frac{1}{\delta })\right).$$

$$\mathcal{O}\left(\frac{\Vert \mathbf{H}{\Vert }_F\Vert {\mathbf{B}}_t\Vert \Vert {\tilde{\mathbf{v}}}^{(i)}\Vert }{{ϵ}^2}(L(\mathbf{H})+Q({\mathbf{B}}_t)+Q(\tilde{\mathbf{V}}))\log (\frac{1}{\delta })\right).$$

$$\mathcal{O}\left(\frac{k\Vert {\hat{\mathbf{q}}}_t{\Vert }^2}{\Vert \tilde{\mathbf{V}}{\hat{\mathbf{q}}}_t{\Vert }^2}(S(\tilde{\mathbf{V}})+kQ(\tilde{\mathbf{V}})\right).$$

$${\mathcal{A}}_t=\arg \underset{i\in [n_{\text{scalebox}0.49Y}]}{\min }{\left({\mathcal{P}}_{{\hat{\mathbf{Y}}}_t}(i)-{\hat{\xi }}_t\underset{j\in [n_{\text{scalebox}0.49X}]}{\max }({\mathcal{P}}_{{\hat{\mathbf{X}}}_t}(j))\right)}_{+}.$$

$$I_{\text{scalebox}0.49X,1}=\arg \underset{z}{\max }{N}_{\text{scalebox}0.49X,z},I_{\text{scalebox}0.49Y,1}=\arg \underset{z}{\max }{N}_{\text{scalebox}0.49Y,z},$$

$$J_{\text{scalebox}0.49X,1}=\arg \underset{z}{\max }{\mathcal{P}}_{{\hat{\mathbf{X}}}_t}(\mathbf{x}=z),\text{and}{J}_{\text{scalebox}0.49Y,1}=\arg \underset{z}{\max }{\mathcal{P}}_{{\hat{\mathbf{Y}}}_t}(\mathbf{y}=z).$$

$${\mathcal{P}}_{{\hat{\mathbf{X}}}_t}(\mathbf{x}=I_{\text{scalebox}0.49X,1})=\underset{J_{X,v}\in [w_{\text{scalebox}0.49X}]}{\max }({\mathcal{P}}_{{\hat{\mathbf{X}}}_t}(J_{\text{scalebox}0.49X,v})),\frac{N_{\text{scalebox}0.49Y,I_{\text{scalebox}0.49Y,v^{\ast }}}}{N_{\text{scalebox}0.49Y}}={\mathcal{P}}_{{\hat{\mathbf{Y}}}_t}(\mathbf{y}={\mathcal{A}}_t),w_{\text{scalebox}0.49X}=n_{\text{scalebox}0.49X},w_{\text{scalebox}0.49Y}=n_{\text{scalebox}0.49Y}.$$

$${\hat{\xi }}_t=\frac{\Vert {\hat{\mathbf{Y}}}_t\Vert }{\Vert {\hat{\mathbf{X}}}_t\Vert }\le {\kappa }_{\text{scalebox}0.49X}\Vert \hat{\mathbf{Y}}{\Vert }_F.$$

$$\mathcal{P}\left(\sum _{i=1}^l\left|\frac{N_i}{N}-p_i\right|\ge \lambda \right)\le 2^l\exp \left(\frac{-N{\lambda }^2}{2}\right).$$

$$\mathcal{P}\left(\left|\frac{N_1}{N_{\text{scalebox}0.49X}}-{\mathcal{P}}_{{\hat{\mathbf{X}}}_t}({\mathcal{X}}_1)\right|+\left|\left(\frac{N_2}{N_{\text{scalebox}0.49X}}-{\mathcal{P}}_{{\hat{\mathbf{X}}}_t}({\mathcal{X}}_2)\right)\right|\ge \lambda \right)\le 4\exp \left(\frac{-N_{\text{scalebox}0.49X}{\lambda }^2}{2}\right).$$

$$\left|\frac{N_{\text{scalebox}0.49X,I_{\text{scalebox}0.49X,1}}}{N_{\text{scalebox}0.49X}}-{\mathcal{P}}_{{\hat{\mathbf{X}}}_t}(I_{\text{scalebox}0.49X,1})\right|$$

$$\le$$

$$\mathcal{P}\left(\left|\frac{N_1}{N_{\text{scalebox}0.49Y}}-{\mathcal{P}}_{{\hat{\mathbf{Y}}}_t}({\mathcal{Y}}_1)\right|+\left|\left(\frac{N_2}{N_{\text{scalebox}0.49Y}}-{\mathcal{P}}_{{\hat{\mathbf{Y}}}_t}({\mathcal{Y}}_2)\right)\right|+\left|\left(\frac{N_3}{N_{\text{scalebox}0.49Y}}-{\mathcal{P}}_{{\hat{\mathbf{Y}}}_t}({\mathcal{Y}}_3)\right)\right|\ge \lambda \right)\le 8\exp \left(\frac{-N_{\text{scalebox}0.49Y}{\lambda }^2}{2}\right).$$

$$\left|\frac{N_{\text{scalebox}0.49Y,I_{\text{scalebox}0.49Y,v^{\ast }}}}{N_{\text{scalebox}0.49Y}}-{\mathcal{P}}_{{\hat{\mathbf{Y}}}_t}(I_{\text{scalebox}0.49Y,v^{\ast }})\right|$$

$$\le$$

$$\Vert {\mathcal{P}}_{{\hat{\mathbf{H}}}_t}-{\mathcal{P}}_{{\mathbf{H}}_t}{\Vert }_{TV}\le \Vert {\mathcal{P}}_{{\hat{\mathbf{H}}}_t}-{\mathcal{P}}_{{\tilde{\mathbf{H}}}_t}{\Vert }_{TV}+\Vert {\mathcal{P}}_{{\tilde{\mathbf{H}}}_t}-{\mathcal{P}}_{{\mathbf{H}}_t}{\Vert }_{TV},$$

$$\Vert {\mathcal{P}}_{{\tilde{\mathbf{H}}}_t}-{\mathcal{P}}_{{\mathbf{H}}_t}{\Vert }_{TV}\le \frac{ϵ}{2}.$$

$$\Vert {\tilde{\mathbf{H}}}_t-{\mathbf{H}}_t\Vert \le \frac{ϵ}{4}\Vert {\mathbf{H}}_t\Vert .$$

Full text

1 Introduction

Maximum a posteriori (MAP) estimation is a central problem in machine and statistical learning [5, 22]. The general MAP problem has been proven to be NP hard [33]. Despite the hardness in the general case, there are two fundamental learning models, the matrix factorization and the latent variable model, that enable MAP problem to be solved in polynomial runtime under certain constraints [24, 25, 29, 30, 32]. The algorithms that have been developed for these learning models have been used extensively in machine learning with competitive performance, particularly on tasks such as subspace clustering, topic modeling, collaborative filtering, structure prediction, feature engineering, motion segmentation, sequential data analysis, and recommender systems [15, 25, 30]. A recent study demonstrates that MAP problems addressed by matrix factorization and the latent variable models can be reduced to the general minimum conical hull problem [41]. In particular, the general minimum conical hull problem transforms problems resolved by these two learning models into a geometric problem, whose goal is to identify a set of extreme data points with the smallest cardinality in dataset $\bm{Y}$ such that every data point in dataset $\bm{X}$ can be expressed as a conical combination of the identified extreme data points. Unlike the matrix factorization and the latent variable models that their optimizations generally suffer from the local minima, a unique global solution is guaranteed for the general minimum conical hull problem [41]. Driven by the promise of a global solution and the broad applicability, it is imperative to seek algorithms that can efficiently resolve the general minimum conical hull problem with theoretical guarantees.

The divide-and-conquer anchoring (DCA) scheme is among the best currently known solutions for addressing general minimum conical hull problems. The idea is to identify all $k$ extreme rays (i.e., $k$ extreme data points) of a conical hull from a finite set of real data points with high probability [41]. The discovered extreme rays form the global solution for the problem with explainability and, thus, the scheme generalizes better than conventional algorithms, such as expectation-maximization [10]. DCA’s strategy is to decompose the original problem into distinct subproblems. Specifically, the original conical hull problem is randomly projected on different low-dimensional hyperplanes to ease computation. Such a decomposition is guaranteed by the fact that the geometry of the original conical hull is partially preserved after a random projection. However, a weakness of DCA is that it requires a polynomial runtime complexity with respect to the size of the input. This complexity heavily limits DCA’s use in many practical situations given the number of massive-scale datasets that are now ubiquitous [38]. Hence, more effective methods for solving general minimum conical hull problems are highly desired.

To address the above issue, we propose an efficient classical algorithm that tackles general minimum conical hull problems in polylogarithmic time with respect to the input size. Consider two datasets $\bm{X}$ and $\bm{Y}$ that have stored in a specific low-overhead data structure, i.e., a sampled-based data structure supports the length-square sampling operations [34]. Let the maximum rank, the Frobenius norm, and the condition number of the given two datasets be $k$, $\|\bm{H}\|_{F}$, and $\kappa$, respectively. We prove that the runtime complexity of the our algorithm is $\tilde{\mathcal{O}}(k^{6}\kappa^{12}\|\bm{H}\|_{F}^{6}/\epsilon^{6})$ with the tolerable level of error $\epsilon$. The achieved sublinear runtime complexity indicates that our algorithm has capability to benefit numerous learning tasks that can be mapped to the general minimum conical hull problem, e.g., the MAP problems addressed by latent variable models and matrix factorization.

Two core ingredients of the proposed algorithm are the ‘divide-and-conquer’ strategy and the reformulation of the minimum conical hull problem as a sampling problem. We adopt the ‘divide-and-conquer’ strategy to acquire a favorable property from DCA. In particular, all subproblems, i.e., the general minimum conical hull problems on different low-dimensional random hyperplanes, are independent of each other. Therefore, they can be processed in parallel. In addition, the total number of subproblems is only polynomially proportional to the rank of the given dataset [40, 41]. An immediate observation is that the efficiency of solving each subproblem governs the efficiency of tackling the general minimum conical hull problem. To this end, our algorithm converts each subproblem into an approximated sampling problem and obtains the solution in sublinear runtime complexity. Through advanced sampling techniques [16, 34], the runtime complexity to prepare the approximated sampling distribution that corresponds to each subproblem is polylogarithmic in the size of input. To enable our algorithm has an end-to-end sublinear runtime, we propose a general heuristic post-selection method to efficiently sample the solution from the approximated distribution, whose computation cost is also polylogarithmic in size.

Our work creates an intriguing aftermath for the quantum machine learning community. The overarching goal of quantum machine learning is to develop quantum algorithms that quadratically or even exponentially reduce the runtime complexity of classical algorithms [4]. Numerous quantum machine learning algorithms with provably quantum speedups have been proposed in the past decade [17, 20, 28]. However, the polylogarithmic runtime complexity in our algorithm implies that a rich class of quantum machine-learning algorithms do not, in fact, achieve these quantum speedups. More specifically, if a quantum algorithm aims to solve a learning task that can be mapped to the general minimum conical hull problem, its quantum speedup will collapse. In our examples, we show that quantum speedups collapse for these quantum algorithms: recommendation system [21], matrix factorization [13], and clustering [1, 27, 37].

Related Work. We make the following comparisons with previous studies in maximum a posteriori estimation and quantum machine learning.

The mainstream methods to tackle MAP problems prior to the study [41] can be separated into two groups. The first group includes expectation-maximization, sampling methods, and matrix factorization [9, 14, 31], where the learning procedure has severely suffered from local optima. The second group contains the method of moments [3], which has suffered from the large variance and may lead to the failure of final estimation. The study [41] effectively alleviates the difficulties encountered by the above two groups. However, a common weakness possessed by all existing methods is the polynomial runtime with respect to the input size. 2. 2.

There are several studies that collapse the quantum speedups by proposing quantum-inspired classical algorithms [7, 8, 16]. For example, the study [34] removes the quantum speedup for recommendation systems tasks; the study [35] eliminates the quantum speedup for principal component analysis problems; the study [11] collapses the quantum speedup for support vector machine problem. The correctness of a branch of quantum-inspired algorithms is validated by study [2]. The studies [34] and [11] can be treated as a special case of our result, since both recommendation systems tasks and support vector machine can be efficiently reduced to the general conical hull problems [26, 40]. In other words, our work is a more general methodology to collapse the speedups achieved by quantum machine learning algorithms.

The rest of this paper proceeds as follows. A formal outline of the general minimum conical hull problem is given in Section 2. The computational complexity of our algorithm is explained and analyzed in Section 3. Section 4 discusses the algorithm’s correctness. We conclude the paper in Section 5.

2 Problem Setup

2.1 Notations

Throughout this paper, we use the following notations. We denote $\{1,2,...,n\}$ as $[n]$. Given a vector ${\bm{v}}\in\mathbb{R}^{n}$, ${\bm{v}}_{i}$ or ${\bm{v}}(i)$ represents the $i$-th entry of $\bm{v}$ with $i\in[n]$ and $\|{\bm{v}}\|$ refers to the $\ell_{2}$ norm of $\bm{v}$ with $\|{\bm{v}}\|=\sqrt{\sum_{i=1}^{n}{\bm{v}}_{i}^{2}}$. The notation $\bm{e}_{i}$ always refers to the $i$-th unit basis vector. Suppose that ${\bm{v}}$ is nonzero, we define $\mathcal{P}_{{\bm{v}}}$ as a probability distribution in which the index $i$ of ${\bm{v}}$ will be chosen with the probability $\mathcal{P}_{{\bm{v}}}(i)=(|{\bm{v}}_{i}|/\|{\bm{v}}\|)^{2}$ for any $i\in[n]$. A sample from $\mathcal{P}_{{\bm{v}}}$ refers to an index number $i$ of $\bm{v}$, which will be sampled with the probability $\mathcal{P}_{{\bm{v}}}(i)$. Given a matrix $\bm{X}\in\mathbb{R}^{n\times m}$, $\bm{X}_{ij}$ represents the $(i,j)$-entry of $\bm{X}$ with $i\in[n]$ and $j\in[m]$. $\bm{X}(i,:)$ and $\bm{X}(:,j)$ represent the $i$-th row and the $j$-th column of the matrix $\bm{X}$, respectively. The transpose of a matrix $\bm{X}$ (a vector ${\bm{v}}$) is denoted by $\bm{X}^{\top}$ (${\bm{v}}^{\top}$). The Frobenius and spectral norm of $\bm{X}$ is denoted as $\|\bm{X}\|_{F}$ and $\|\bm{X}\|_{2}$, respectively.The condition number $\kappa_{\scalebox{0.49}{X}}$ of a positive semidefinite $\bm{X}$ is $\kappa_{\scalebox{0.49}{X}}=\|\bm{X}\|_{2}/\sigma_{\min}(\bm{X})$, where $\sigma_{\min}(\bm{X})$ refers to the minimum singular of $\bm{X}$. $\mathbb{R}_{+}$ refers to the real positive numbers. Given two sets $\mathcal{A}$ and $\mathcal{B}$, we denote $\mathcal{A}$ minus $\mathcal{B}$ as $\mathcal{A}\setminus\mathcal{B}$. The cardinality of a set $\mathcal{A}$ is denoted as $|\mathcal{A}|$. We use the notation $\tilde{\mathcal{O}}(k)$ as shorthand for $\mathcal{O}(k\log(n))$.

2.2 General minimum conical hull problem

A cone is a non-empty convex set that is closed under the conical combination of its elements. Mathematically, given a set $\mathcal{R}=\{r_{i}\}_{i=1}^{k}$ with $r_{i}$ being the ray, a cone formed by $\mathcal{R}$ is defined as ${cone}(\mathcal{R})=\{\sum_{i=1}^{k}\alpha_{i}r_{i}:~{}\forall i\in[k],~{}r_{i}\in\mathcal{R},~{}\alpha_{i}\in\mathbb{R}_{+}\}$, and ${cone}(\mathcal{R})$ is the conical hull of $\mathcal{R}$. A ray $r_{i}$ is an extreme ray (or an anchor) if it cannot be expressed as the conical combination of elements in $\mathcal{R}\setminus r_{i}$. A fundamental property of the cone and conical hull is its separability, namely, whether a point in ${cone}(\mathcal{R})$ can be represented as a conical combination of certain subsets of rays that define ${cone}(\mathcal{R})$. The above separability can be generalized to the matrix form as follows [41].

Definition 1** (General separability condition [41]).**

Let $\bm{X}\in\mathbb{R}^{n_{\scalebox{0.49}{X}}\times m}$ be a matrix with $n_{\scalebox{0.49}{X}}$ rows. Let $\bm{Y}\in\mathbb{R}^{n_{\scalebox{0.49}{Y}}\times m}$ be a matrix with $n_{\scalebox{0.49}{Y}}$ rows, and let $\bm{Y}_{\mathcal{A}}$ be a submatrix of $\bm{Y}$ with rows ${\mathcal{A}}\subset[n_{\scalebox{0.49}{Y}}]$. We say that $\bm{X}$ is separable with respect to $\bm{Y}_{\mathcal{A}}$ if $\bm{X}=\bm{F}\bm{Y}_{\mathcal{A}}$, where $rank(\bm{F})\geq|\mathcal{A}|$.

In other words, the general separability condition states that, $\forall i\in[n_{\scalebox{0.49}{X}}]$, we have $\bm{X}_{i}\in{cone}(\bm{Y}_{\mathcal{A}})$, where the set $\bm{Y}_{\mathcal{A}}=\{\bm{Y}(i;:)\}_{i\in{\mathcal{A}}}.$ Under this definition, the general minimum conical hull problem aims to find the minimal set ${\mathcal{A}}$, as the so-called anchor set, from the rows of $\bm{Y}$.

Definition 2** (General minimum conical hull problem [41]).**

Given two matrices $\bm{X}\in\mathbb{R}^{n_{\scalebox{0.49}{X}}\times m}$ and $\bm{Y}\in\mathbb{R}^{n_{\scalebox{0.49}{Y}}\times m}$ with $n_{X}$ and $n_{Y}$ rows, respectively, the general minimum conical hull problem ${\rm{MCH}}(\bm{Y},\bm{X})$ finds a minimal subset of rows in $\bm{Y}$ whose conical hull contains $cone(\bm{X})$:

[TABLE]

where $\text{cone}(\bm{Y}_{\mathcal{A}})$ is induced by rows ${\mathcal{A}}$ of $\bm{Y}$.

We remark that the general separability condition is reasonable for many learning tasks, i.e., any learning task can be solved by the matrix factorization model or the latent variable model possessing general separability [41].

2.3 Divide-and-conquer anchoring scheme for the general minimum conical hull problem

Efficiently locating the anchor set ${\mathcal{A}}$ of a general minimum conical hull problem is challenging. To resolve this issue, the divide-and-conquer anchoring learning (DCA) scheme is proposed [40, 41]. This scheme adopts the following strategy. The original minimum conical hull problem is first reduced to a small number of subproblems by projecting the original cone into low-dimensional hyperplanes. Then these subproblems on low-dimensional hyperplanes are then tackled in parallel. The anchor set of the original problem can be obtained by combining the anchor sets of the subproblems because the geometric information of the original conical hull is partially preserved after projection. Moreover, the efficiency of DCA schemes can be guaranteed because anchors in these subproblems can be effectively determined in parallel and the total number of subproblems is modest.

DCA is composed of two steps, i.e., the divide step and the conquer step. Following the notations of Definition 2, given two sets of points (rows) with $\bm{X}\in\mathbb{R}^{n_{\scalebox{0.49}{X}}\times m}$ and $\bm{Y}\in\mathbb{R}^{n_{\scalebox{0.49}{Y}}\times m}$, the goal of DCA is to output an anchor set ${\mathcal{A}}$ such that $\bm{X}(i;:)\in cone(\bm{Y}_{\mathcal{A}})$ with ${\mathcal{A}}\subset[n_{\scalebox{0.49}{Y}}]$ and $|{\mathcal{A}}|=k$, $\forall i\in[n_{\scalebox{0.49}{X}}]$.

Divide step. A set of projection matrices $\{\bm{B}_{t}\}_{t=1}^{p}$ with $\bm{B}_{t}\in\mathbb{R}^{m\times d}$ is sampled from a random ensemble $\mathbb{M}$, e.g., Gaussian random matrix ensemble, the ensemble composed of standard unit vectors $\bm{e}_{i}$ or real data vectors, and various sparse random matrix ensembles [41]. The general minimum conical hull problem for the $t$-th subproblem is to find an anchor set ${\mathcal{A}}_{t}$ for the projected matrices $\bm{Y}_{t}:=\bm{Y}\bm{B}_{t}\in\mathbb{R}^{n_{\scalebox{0.49}{Y}}\times d}$ and $\bm{X}_{t}:=\bm{X}\bm{B}_{t}\in\mathbb{R}^{n_{\scalebox{0.49}{X}}\times d}$:

[TABLE]

where ${\bm{Y}}_{{\mathcal{A}}_{t}}$ is the submatrix of $\bm{Y}$ whose rows are indexed by ${\mathcal{A}}_{t}\subset[n_{\scalebox{0.49}{Y}}]$.

Conquer step. This step yields the anchor set ${\mathcal{A}}$ by employing the following selection rule to manipulate the collected $\{{{\mathcal{A}}}_{t}\}_{t=1}^{p}$. First, we compute, $\forall i\in[n_{\scalebox{0.49}{Y}}]$,

[TABLE]

where $\bm{Y}(i;:)$ is the $i$-th row of $\bm{Y}$, and $\mathbbm{1}_{{{\mathcal{A}}}_{t}}(\bm{Y}(i;:))$ is the indicator function that outputs ‘1’ when the index $i$ is in ${{\mathcal{A}}}_{t}$, and zero otherwise. The anchor set ${\mathcal{A}}$ of size $k$ is constructed by selecting $k$ indexes with the largest $\hat{g}_{i}$.

It has been proved, with high probability that, solving $p=\mathcal{O}(k\log k)$ subproblems are sufficient to find all anchors in $\mathcal{A}$, where $k$ refers to the number of anchors [41]. The total runtime complexity of DCA is $\tilde{\mathcal{O}}(\max\{\text{poly}(n_{\scalebox{0.49}{X}}),\text{poly}(n_{\scalebox{0.49}{Y}})\})$ when parallel processing is allowed.

2.4 Reformulation as a sampling problem

Interestingly, each subproblem in Eqn. (2) can be reduced to a sampling problem when $d=1$. This observation is one of the crucial components that make our algorithm exponentially faster than the conventional DCA schemes [41].

Fix $d=1$. Following the result of [39], Eqn. (2) becomes

[TABLE]

where $\bm{Y}_{t}\in\mathbb{R}^{n_{\scalebox{0.49}{Y}}\times 1}$, $\bm{X}_{t}\in\mathbb{R}^{n_{\scalebox{0.49}{X}}\times 1}$, $(x)_{+}=x$ if $x\geq 0$ and $\infty$ otherwise. We give an intuitive explanation of Eqn. (4) in Figure (1).

Note that Eqn. (4) can then be written as

[TABLE]

where $\mathcal{P}_{\bm{X}_{t}}$ and $\mathcal{P}_{\bm{Y}_{t}}$ refer to the distributions of $\bm{X}_{t}$ and $\bm{Y}_{t}$, and $\xi_{t}=\|\bm{X}_{t}\|/\|\bm{Y}_{t}\|$ is a constant at $t$-th random projection to rescale the value $\max_{j\in[n_{\scalebox{0.49}{X}}]}(\mathcal{P}_{\bm{X}_{t}}(j))$ . We will explain how to efficiently approximate $\xi_{t}$ in Section 3.2.

3 Main Algorithm

Our algorithm generates two distributions $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ to approximate the targeted distributions $\mathcal{P}_{\bm{X}_{t}}$ and $\mathcal{P}_{\bm{Y}_{t}}$ as defined in Eqn. (5) for any $t\in[p]$. The core of our algorithm consists of the following major steps. In the preprocessing step, we reconstruct two matrices, $\tilde{\bm{X}}$ and $\tilde{\bm{Y}}$, to approximate the original matrices $\bm{X}$ and $\bm{Y}$ so that the $t$-th projected subproblem can be replaced with $\tilde{\bm{X}}_{t}=\tilde{\bm{X}}\bm{B}_{t}$ and $\bm{\tilde{Y}}_{t}=\tilde{\bm{Y}}\bm{B}_{t}$ with little disturbance. The main tool for this approximation is the subsampling method with the support of the square-length sampling operations [34]. This step also appears in [8, 16, 34]. All subproblems are processed in parallel, following the divide-and-conquer principle. The divide step employs two sampling subroutines that allow us to efficiently generate two distributions $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ to approximate $\mathcal{P}_{\bm{\tilde{X}}_{t}}$ and $\mathcal{P}_{\tilde{\bm{Y}}_{t}}$ (equivalently, $\mathcal{P}_{\bm{{X}}_{t}}$ and $\mathcal{P}_{{\bm{Y}}_{t}}$). We then propose the general heuristic post-selection rule to identify the target index $\mathcal{A}_{t}$ by substituting $(\mathcal{P}_{\bm{Y}_{t}}(i)-\xi_{t}\max_{j\in[n_{\scalebox{0.49}{X}}]}(\mathcal{P}_{\bm{X}_{t}}(j))_{+}$ in Eqn (5) with $(\mathcal{P}_{\bm{\hat{Y}}_{t}}(i)-\xi_{t}\max_{j\in[n_{\scalebox{0.49}{X}}]}(\mathcal{P}_{\bm{\hat{X}}_{t}}(j))_{+}$. Lastly, we employ the selection rule in Eqn. (3) to form the anchor set $\mathcal{A}$ for the original minimum conical hull problem.

Before elaborating the details of the main algorithm, we first emphasize the innovations of this work, i.e., the reformulation of the general conical hull problem as a sampling problem and the general heuristic post-selection rule. The sampling version of the general conical hull problem is the precondition to introduce advanced sampling techniques to reduce the computational complexity. The general heuristic post-selection rule is the central component to guarantee that the solution of the general conical hull problem can be obtained in sublinear runtime. In particular, the intrinsic mechanism of the general conical hull problem enables us to employ the general heuristic post-selection rule to query a specific element from the output in polylogarithmic runtime.

In this section, we introduce the length-square sampling operations in Subsection 3.1. The implementation of our algorithm is shown in Subsection 3.2. The computation analysis is given in Subsection 3.3. We give the proof of Theorems 5 and 6 in supplementary material C and E, respectively.

3.1 Length-square sampling operations

The promising performance of various sampling algorithms for machine learning tasks is guaranteed when the given dataset supports length-square sampling operations [8, 16, 34]. We first give the definition of the access cost to facilitate the description of such sampling operations,

Definition 3** (Access cost).**

Given a matrix $\bm{W}\in\mathbb{R}^{n\times m}$, we denote that the cost of querying an entry $\bm{W}(i,j)$ or querying the Frobenius norm $\|\bm{W}\|_{F}$ is $Q(\bm{W})$, the cost of querying the $\ell_{2}$ norm of $\|\bm{W}(i,:)\|$ is $N(\bm{W})$, and the cost of sampling a row index $i\in[n]$ of $\bm{H}$ with the probability $\mathcal{P}_{\bm{W}}(i)=\|\bm{W}(i,:)\|^{2}/\|\bm{W}\|_{F}^{2}$ or sampling an index $j\in[m]$ with the probability $\mathcal{P}_{\bm{W}(i,:)}(j)$ is $S(\bm{W})$. We denote the overall access cost of $\bm{W}$ as $L(\bm{W}):=S(\bm{W})+Q(\bm{W})+N(\bm{W})$.

We use an example to address the difference between query access and sampling access. Given a vector $\bm{v}\in\mathbb{R}^{n}$, the problem is to find a hidden large entry $\bm{v}(i)$. It takes $\Omega(n)$ cost to find $\bm{v}(i)$ with just query access, while the computation cost is $\mathcal{O}(1)$ with query and sample access.

The length-square sampling operations, as so-called $\ell_{2}$ norm sampling operations, are defined as:

Proposition 4** ($\ell_{2}$ norm sampling operation).**

Given an input matrix $\bm{H}\in\mathbb{R}^{n\times m}$ with $s$ non-zero entries, there exists a data structure storing $\bm{H}$ in space $\mathcal{O}(s\log^{2}m)$, with the following properties:

The query cost $Q(\bm{H})$ is at most $\mathcal{O}(\log(nm))$; 2. 2.

The query cost $N(\bm{H})$ is at most $\mathcal{O}(\log^{2}(n))$; 3. 3.

The sampling cost $S(\bm{H})$ is at most $\mathcal{O}(\log^{2}(nm))$.

The overall cost of accessing $\bm{H}$ is therefore $L(\bm{H})=\mathcal{O}(poly(\log(mn)))$.

Remark. The $\ell_{2}$ norm sampling operation can be efficiently fulfilled if the input data are stored in a low-overhead data structure, e.g., the binary tree structure (BNS) [21] (more details about BNS are given in supplementary material A).

3.2 The implementation of the algorithm

Our algorithm consists of three steps, the preprocessing step, the divide step, and the conquer step. The first step prepares an efficient description for $\bm{X}$ and $\bm{Y}$. The second step locates anchors $\{\mathcal{A}_{t}\}_{t=1}^{p}$ by solving $p$ subproblems in parallel. The last step obtains the anchor set $\mathcal{A}$.

Preprocessing step. The preprocessing step aims to efficiently construct $\bm{\tilde{X}}$ and $\bm{\tilde{Y}}$ such that the matrix norm $\|\bm{\tilde{X}}-\bm{X}\|_{2}$ and $\|\bm{\tilde{Y}}-\bm{Y}\|_{2}$ are small. This step employs the subsampling method [34] to construct an approximated left singular matrix $\tilde{\bm{V}}_{\scalebox{0.49}{H}}$ of $\bm{H}$, where $\bm{H}\in\mathbb{R}^{n_{\scalebox{0.49}{H}}\times m_{\scalebox{0.49}{H}}}$ can be either $\bm{X}$ or $\bm{Y}$, so that $\bm{\tilde{H}}=\tilde{\bm{V}}_{\scalebox{0.49}{H}}\tilde{{\bm{V}}_{\scalebox{0.49}{H}}}^{\top}\bm{H}$. If no confusion arises, the subscript $\bm{H}$ can be disregarded.

We summarize the subsampling method in Algorithm 1 to detail the acquisition of $\tilde{\bm{V}}$. We build the matrix $\bm{R}\in\mathbb{R}^{n\times s}$ by sampling $s$ columns from $\bm{H}$ and then build the matrix $\bm{C}\in\mathbb{R}^{s\times s}$ by sampling $s$ rows from $\bm{R}$. After obtaining $\bm{R}$ and $\bm{C}$, we implicitly define the approximated left singular matrix $\tilde{\bm{V}}\in\mathbb{R}^{n\times k}$ as

[TABLE]

where $\{\sigma^{(i)}\}_{i=1}^{k}$ and $\{\bm{\omega}^{(i)}\}_{i=1}^{k}$ refer to $k$ singular values and right singular vectors of $\bm{C}$111The ‘implicitly define $\bm{\tilde{V}}$’ means that only the index array of $\bm{R}$ and $\bm{C}$, and the SVD result of $\bm{C}$ are required to be stored in the memory..

Divide step. The obtained approximated left singular matrices $\tilde{\bm{V}}_{\scalebox{0.49}{X}}\in\mathbb{R}^{n_{\scalebox{0.49}{X}}\times k_{\scalebox{0.49}{X}}}$ and $\tilde{\bm{V}}_{\scalebox{0.49}{Y}}\in\mathbb{R}^{n_{\scalebox{0.49}{Y}}\times k_{\scalebox{0.49}{Y}}}$ enable us to employ advanced sampling techniques to locate potential anchors $\{\mathcal{A}_{t}\}_{t=1}^{p}$. Here we only focus on locating the anchor $\mathcal{A}_{t}$ for the $t$-th subproblem, since each subproblem is independent and can be solved in the same way. The divide step employs Eqn. (5) to locate $\mathcal{A}_{t}$. In particular, we first prepare two distributions $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ to approximate $\mathcal{P}_{\bm{X}_{t}}$ and $\mathcal{P}_{\bm{Y}_{t}}$, and then sample these two distributions to locate $\mathcal{A}_{t}$.

The preparation of two distributions $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ is achieved by exploiting two sampling subroutines, the inner product subroutine and the thin matrix-vector multiplication subroutine [34]. We detail these two subroutines in supplementary material B.2. Recall that the two approximated matrices at the $t$-th subproblem are ${\bm{\tilde{X}}_{t}}=\tilde{{\bm{V}}}_{\scalebox{0.49}{X}}(\tilde{{\bm{V}}}_{\scalebox{0.49}{X}}^{\top}\bm{X}_{t})$ and ${\bm{\tilde{Y}}_{t}}=\tilde{\bm{V}}_{\scalebox{0.49}{Y}}(\tilde{{\bm{V}}}_{\scalebox{0.49}{Y}}^{\top}\bm{Y}_{t}).$ Denote $\bm{\tilde{q}}_{\scalebox{0.49}{X},t}\equiv\tilde{{\bm{V}}}_{\scalebox{0.49}{X}}^{\top}\bm{X}_{t}\in\mathbb{R}^{k_{\scalebox{0.49}{X}}\times 1}$ and $\bm{\tilde{q}}_{\scalebox{0.49}{Y},t}\equiv\tilde{{\bm{V}}}_{\scalebox{0.49}{Y}}^{\top}\bm{Y}_{t}\in\mathbb{R}^{k_{\scalebox{0.49}{Y}}\times 1}$. Instead of directly computing $\bm{\tilde{q}}_{\scalebox{0.49}{X},t}$ and $\bm{\tilde{q}}_{\scalebox{0.49}{Y},t}$, we construct their approximated vectors $\hat{\bm{q}}_{\scalebox{0.49}{X},t}$ and $\hat{\bm{q}}_{\scalebox{0.49}{Y},t}$ using the inner product subroutine to ensure the low computational cost, followed by the thin matrix-vector multiplication subroutine to prepare probability distributions $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$, where $\bm{\hat{X}}_{t}\equiv\tilde{{\bm{V}}}_{\scalebox{0.49}{X}}\hat{\bm{q}}_{\scalebox{0.49}{X},t}$ and $\bm{\hat{Y}}_{t}\equiv\tilde{{\bm{V}}}_{\scalebox{0.49}{Y}}\hat{\bm{q}}_{\scalebox{0.49}{Y},t}$. The closeness between $\mathcal{P}_{\bm{\hat{X}}_{t}}$ (resp. $\mathcal{P}_{\bm{\hat{Y}}_{t}}$) and $\mathcal{P}_{\bm{{X}}_{t}}$ (resp. $\mathcal{P}_{\bm{{Y}}_{t}}$) is controlled by the number of samplings $s$, as analyzed in Section 4.

The rescale parameter $\xi_{t}$ defined in Eqn. (5) can be efficiently approximated by employing the inner product subroutine. Recall that $\xi_{t}=\|\bm{Y}_{t}\|/\|\bm{X}_{t}\|$. We can approximate $\xi_{t}$ by $\hat{\xi}_{t}=\|\bm{\hat{Y}}_{t}\|/\|\bm{\hat{X}}_{t}\|$. Alternatively, an efficient method of approximating $\|\bm{\hat{Y}}_{t}\|$ and $\|\bm{\hat{X}}_{t}\|$ is sufficient to acquire $\xi_{t}$. Let $\bm{H}$ be the general setting that can either be $\bm{X}$ or $\bm{Y}$. The $\ell_{2}$ norm of $\|\bm{\hat{H}}_{t}\|$ can be efficiently estimated by using the inner product subroutine. Intuitively, the $\ell_{2}$ norm of $\|\hat{\bm{H}}_{t}\|$ can be expressed by the inner product of $\hat{\bm{H}}_{t}$, i..e, $\|\hat{\bm{H}}_{t}\|^{2}={\hat{\bm{H}}_{t}^{\top}\hat{\bm{H}}_{t}}$. Recall that $\hat{\bm{H}}_{t}$ has an explicit representation $\hat{\bm{H}}_{t}=\tilde{\bm{V}}_{\scalebox{0.49}{H}}\hat{\bm{q}}_{\scalebox{0.49}{H,t}}$, and the efficient access cost for $\tilde{\bm{V}}_{\scalebox{0.49}{H}}$ and $\tilde{\bm{q}}_{\scalebox{0.49}{H,t}}$ enables us to use the inner product subroutine to efficiently obtain $\|\bm{\tilde{H}}_{t}\|$.

Given $\mathcal{P}_{\bm{\hat{X}}_{t}}$, $\mathcal{P}_{\bm{\hat{Y}}_{t}}$, and $\hat{\xi_{t}}$, we propose the general heuristic post-selection method to determine $\mathcal{A}_{t}$. Following the sampling version of the general minimum conical hull problem defined in Eqn. (5), we first sample the distribution $\mathcal{P}_{\bm{\hat{X}}_{t}}$ with $N_{\scalebox{0.49}{X}}$ times to obtain a value $C_{\bm{\hat{X}}_{t}}^{*}$ such that $C_{\bm{\hat{X}}_{t}}^{*}=\max_{j\in[n_{\scalebox{0.49}{X}}]}(\mathcal{P}_{\bm{\hat{X}}_{t}}(j))$. We next sample the distribution $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ with $N_{\scalebox{0.49}{Y}}$ times to find an index $\hat{\mathcal{A}}_{t}$ approximating ${\mathcal{A}}_{t}=\arg\min_{i\in[n_{\scalebox{0.49}{Y}}]}(\mathcal{P}_{\bm{\hat{Y}}_{t}}(i)-\hat{\xi}_{t}C_{\bm{\hat{X}}_{t}}^{*})_{+}$. The following theorem quantifies the required number of samplings to guarantee $\hat{\mathcal{A}}_{t}=\mathcal{A}_{t}$, where $\mathcal{A}_{t}$ is defined in Eqn. (5).

Theorem 5** (General heuristic post-selection (Informal)).**

Assume that $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ are multinomial distributions. Denote that $C_{\bm{\hat{X}}_{t}}^{*}\geq\varepsilon_{T}$. If $|\mathcal{P}_{\bm{\hat{X}}_{t}}(i)-\mathcal{P}_{\bm{\hat{X}}_{t}}(j)|>\varepsilon$ for $C_{\bm{\hat{X}}_{t}}^{*}=\mathcal{P}_{\bm{\hat{X}}_{t}}(i)$ and $\forall j\neq i$, and $|\mathcal{P}_{\bm{\hat{Y}}_{t}}(\mathcal{A}_{t})-\mathcal{P}_{\bm{\hat{Y}}_{t}}(j)|>\varepsilon$ for $\forall j\neq\mathcal{A}_{t}$ and a small constant $\varepsilon$, then for any $\delta>0$ with a probability at least $1-\delta$, we have $\hat{\mathcal{A}}_{t}=\mathcal{A}_{t}$ with $N_{\scalebox{0.49}{X}},N_{\scalebox{0.49}{Y}}\sim\mathcal{O}(\kappa_{\scalebox{0.49}{X}},polylog(\max\{n_{\scalebox{0.49}{X}},n_{\scalebox{0.49}{Y}}\}))$.

Conquer step.* * After the divide step, a set of potential anchors $\{\mathcal{A}_{t}\}_{t=1}^{p}$ with $p=\mathcal{O}(k\log(k))$, where $k$ refers to the number of anchors with $k\sim\min\{\log(n_{\scalebox{0.49}{X}}),\log(n_{\scalebox{0.49}{Y}})\}$, is collected [13, 39, 41]. We adopt the selection rule defined in Eqn. (3) to determine the anchor set $\mathcal{A}$. This step can be accomplished by various sorting algorithms with runtime $\mathcal{O}(poly(k\log(k)))$ [41].

We summarize our algorithm as follows.

3.3 Computation complexity of the algorithm

The complexity of the proposed algorithm is dominated by the preprocessing step and the divide step. Specifically, the computational complexity is occupied by four elements: finding the left singular matrix $\tilde{\bm{V}}$ (Line 5 in Algorithm 2), estimating $\tilde{\bm{q}}_{t}$ by $\hat{\bm{q}}_{t}$ (Line 7 in Algorithm 2), preparing the approximated probability distribution $\mathcal{P}_{\hat{\bm{H}}}$ (Line 8 in Algorithm 2), and estimating the rescale factor $\hat{\xi}_{t}$ (Line 10 in Algorithm 2). We evaluate the computation complexity of these four operations separately and obtain the following theorem

Theorem 6**.**

Given two datasets $\bm{X}\in\mathbb{R}^{n_{\scalebox{0.49}{X}}\times m}$ and $\bm{Y}\in\mathbb{R}^{n_{\scalebox{0.49}{Y}}\times m}$ that satisfy the general separability condition and support the length-square sampling operations, the rank and condition number for $\bm{X}$ (resp. $\bm{Y}$) are denoted as $k_{\scalebox{0.49}{X}}$ (resp. $k_{\scalebox{0.49}{Y}}$) and $\kappa_{\scalebox{0.49}{X}}$ (resp. $\kappa_{\scalebox{0.49}{Y}}$). The tolerable error is denoted as $\epsilon$. The runtime complexity of the proposed algorithm for solving the general minimum conical hull problem is

[TABLE]

4 Correctness of The Algorithm

In this section, we present the correctness of our algorithm. We also briefly explain how the results are derived and validate our algorithm by numerical simulations. We provide the detailed proofs of Theorem 7 in the supplementary material D.

The correctness of our algorithm is determined by the closeness between the approximated distribution and the analytic distribution. The closeness is evaluated by the total variation distance [34], i.e., Given two vectors $\mathbf{v}\in\mathbb{R}^{n}$ and $\mathbf{w}\in\mathbb{R}^{n}$, the total variation distance of two distributions $\mathcal{P}_{{\mathbf{v}}}$ and $\mathcal{P}_{{\mathbf{w}}}$ is $\|\mathcal{P}_{{\mathbf{v}}},\mathcal{P}_{{\mathbf{w}}}\|_{TV}:=\frac{1}{2}\sum_{i=1}^{n}|\mathcal{P}_{{\mathbf{v}}}(i)-\mathcal{P}_{{\mathbf{w}}}(i)|$. The following theorem states that $\|\mathcal{P}_{\bm{\hat{X}}_{t}}-\mathcal{P}_{\bm{{X}}_{t}}\|_{TV}$ and $\|\mathcal{P}_{\bm{\hat{Y}}_{t}}-\mathcal{P}_{\bm{{Y}}_{t}}\|_{TV}$ are controlled by the number of samplings $s$:

Theorem 7**.**

Given a matrix $\bm{H}\in\mathbb{R}^{n\times m}$ with $\ell_{2}$ norm sampling operations, let $\bm{R}\in\mathbb{R}^{n\times s}$, $\bm{C}\in\mathbb{R}^{s\times s}$ be the sampled matrices from $\bm{H}$ following Algorithm 1. The distribution prepared by the proposed algorithm is denoted as $\mathcal{P}_{\bm{\hat{H}}_{t}}$. Set

[TABLE]

with probability at least $(1-\eta)$, we always have $\|\mathcal{P}_{\bm{\hat{H}}_{t}}-\mathcal{P}_{\bm{{H}}_{t}}\|_{TV}\leq\epsilon$.

We apply the proposed algorithm to accomplish the near separable nonnegative matrix factorization (SNMF) to validate the correctness of Theorem 7 [40, 23]. SNMF, which has been extensive applied to hyperspectral imaging and text mining, can be treated as a special case of the general minimum conical hull problem [41]. Given a non-negative matrix $\bm{X}$, SNMF aims to find a decomposition such that $\bm{X}=\bm{F}\bm{X}(\mathcal{R},:)+\bm{N}$, where the basis matrix $\bm{X}(\mathcal{R},:)$ is composited by $r$ rows from $\bm{X}$, $\bm{F}$ is the non-negative encoding matrix, and $\bm{N}$ is a noise matrix [12].

The synthetic data matrix used in the experiments is generated in the form of $\bm{Y}\equiv\bm{X}=\bm{F}\bm{X}_{\mathcal{A}}+\bm{N}_{G}$, where $\bm{F}=[\bm{I}_{k};\bm{\mathcal{U}}]$, the entries of noise $\bm{N}_{G}$ are generated by sampling Gaussian distribution $\mathcal{N}(0,\mu)$ with noise level $\mu$, the entries of $\bm{X}_{\mathcal{A}}$ and $\bm{\mathcal{U}}$ are generated by i.i.d. uniform distribution in $[0,1]$ at first and then their rows are normalized to have unit $\ell_{1}$ norm. Analogous to DCA, we use the recovery rate as the metric to evaluate the precision of anchor recovery. The anchor index recovery rate $\rho$ is defined as $\rho=|\mathcal{A}\cap\hat{\mathcal{A}}|/|\mathcal{A}|$, where $\hat{\mathcal{A}}$ refers to the anchor set obtained by our algorithm or DCA.

We set the dimensions of $\bm{X}$ as $500\times 500$ and set $k$ as 10. We generate four datasets with different noise levels, which are $\mu=0,0.1,0.5,2$. The number of subproblems is set as $p=100$. We give nine different settings for the number of sampling $s$ for our algorithm, ranging from $500$ to $8000$. We compare our algorithm with DCA to determine the anchor set $\mathcal{A}$. The simulation results are shown in Figure 2. For the case of $\mu=0$, both our algorithm and DCA can accurately locate all anchors. With increased noise, the recovery rate continuously decreases both for our algorithm and DCA, since the separability assumption is not preserved. As shown in Figure 2 (b), the reconstruction error, which is evaluated by the Frobenius norm of $\|\hat{\bm{X}}-\bm{X}\|_{F}$, continuously decreases with increased $s$ for the noiseless case. In addition, the variance of the reconstructed error, illustrated by the shadow region, continuously shrinks with increased $s$. For the high noise level case, the collapsed separability assumption arises that the reconstruction error is unrelated to $s$.

5 Conclusion

In this paper, we have proposed a sublinear runtime classical algorithm to resolve the general minimum conical hull problem. We first reformulated the general minimum conical hull problem as a sampling problem. We then exploited two sampling subroutines and proposed the general heuristic post-selection method to achieve low computational cost. We theoretically analyzed the correctness and the computation cost of our algorithm. The proposed algorithm benefit numerous learning tasks that can be mapped to the general minimum conical hull problem, especially for tasks that need to process datasets on a huge scale. There are two promising future directions. First, we will explore other advanced sampling techniques to further improve the polynomial dependence in the computation complexity. Second, we will investigate whether there exist other fundamental learning models that can be reduced to the general minimum conical hull problem. One of the most strongest candidates is the semi-definite programming solver.

6 Note Added

Recently we became aware of an independent related work [6], which employs the divide-and-conquer anchoring strategy to solve separable nonnegative matrix factorization problems. Since a major application of the general conical hull problem is to solve matrix factorization, their result can be treated as a special case of our study. Moreover, our study adopts more advanced techniques and provides a better upper complexity bounds than theirs.

Appendix A The Binary Tree Structure for Length-square Sampling Operations

As mentioned in the main text, a feasible solution to fulfill $\ell_{2}$ norm sampling operations is the binary tree structure (BNS) to store data [21]. Here we give the intuition about how BNS constructed for a vector. For ease of notations, we assume the given vector has size $4$ with $\mathbf{v}\in\mathbb{R}^{4}$. As demonstrated in Figure 3, the root node records the square $\ell_{2}$ norm of $\mathbf{v}$. The $i$-th leaf node records the $i$-th entry of $\mathbf{v}(i)$ and its square value, e.g., $\{|\mathbf{v}(i)|^{2},\text{sign}(\mathbf{v}(i))\}$. Each internal node contains the sum of the values of its two immediate children. Such an architecture ensures the $\ell_{2}$ norm sampling opeartions.

Appendix B Two Sampling Subroutines

Here we introduce two sampling subroutines, the inner product subroutine and the thin matrix-vector multiplication subroutine [34], used in the proposed algorithm.

B.1 Inner product subroutine

In our algorithm, the inner product subroutine is employed to obtain each entry of $\hat{\bm{q}}_{t}$ in parallel, i.e., $\hat{\bm{q}}_{t}(i)$ estimates $\tilde{\bm{q}}_{t}(i)\equiv\tilde{\bm{v}}^{(i)\top}\bm{H}_{t}$, where $\tilde{\bm{v}}^{(i)}\in\mathbb{R}^{n\times 1}$, $\bm{H}_{t}\equiv\bm{H}\bm{B}_{t}$, $\bm{H}\in\mathbb{R}^{n\times m}$ and $\bm{B}_{t}\in\mathbb{R}^{m\times 1}$. Let $Z$ be a random variable that, for $j\in[n],l\in[m]$, takes value

[TABLE]

with probability

[TABLE]

We can estimate $\tilde{\bm{q}}_{t}(i)$ using $Z$ as follows [8]. Fix $\epsilon,\delta>0$. Let

[TABLE]

We first sample the distribution $\mathcal{P}_{Z}$ with $N_{Z}$ times to obtain a set of samples $\{z_{i}\}_{i=1}^{N_{Z}}$, followed by dividing them into $N_{p}$ groups, $\{S_{1},\cdots,S_{N_{p}}\}$, where each group $S_{i}$ contains $N_{q}$ samples. Let $\bar{z}_{i}=\sum_{i=1}^{N_{q}}z_{i}/N_{q}$ be the empirical mean of the $i$-th group $S_{i}$, and let $\bar{z}^{*}$ be the median of $\{\bar{z}_{i}\}_{i=1}^{N_{p}}$. Then [16, Lemma 12] and [8, Lemma 7] guarantee that, with probability at least $1-\delta$, the following holds.

[TABLE]

The computational complexity of the inner product subroutine is:

Lemma 8** ([16, Lemma 12] and [8, Lemma 7]).**

Assume that the overall access to $\bm{H}$ is $L(\bm{H})$ and the query access to $\bm{B}_{t}$ and $\tilde{\bm{V}}$ is $Q(\bm{B}_{t})$ and $Q(\tilde{\bm{V}})$, respectively. The runtime complexity to yield Eqn. (9) is

[TABLE]

Proof.

We first recall the main result of Lemma 12 in [16]. Given the overall access to $\bm{H}$ and query access to the matrix $\bm{\tilde{V}}$ and $\bm{B}_{t}$ with complexity $Q(\bm{\tilde{V}})$ and $Q(\bm{B}_{t})$, the inner product $\tilde{\bm{q}}_{t}(i)\equiv\tilde{\bm{v}}^{(i)\top}\bm{H}_{t}$ can be estimated to precision $\epsilon\|\bm{H}\|_{F}\|\tilde{\bm{v}}^{(i)}\|\|\bm{B}_{t}\|_{F}$ with probability at least $1-\delta$ in time

[TABLE]

With setting the precision to $\epsilon$ instead of $\epsilon\|\bm{H}\|_{F}\|\|\tilde{\bm{v}}^{(i)}\|\bm{B}_{t}\|_{F}$, it can be easily inferred that the runtime complexity to estimate $\tilde{\bm{q}}_{t}(i)\equiv\tilde{\bm{v}}^{(i)\top}\bm{H}_{t}$ with probability at least $1-\delta$ is

[TABLE]

∎

B.2 Thin matrix-vector multiplication subroutine

Given a matrix $\tilde{\bm{V}}\in\mathbb{R}^{n\times k}$ and $\bm{\hat{q}}_{t}\in\mathbb{R}^{k}$ with the $\ell_{2}$ norm sampling access, the thin matrix-vector multiplication subroutine aims to output a sample from $\mathcal{P}_{\tilde{\bm{V}}\bm{\hat{q}}_{t}}$. The implementation of the thin matrix-vector multiplication subroutine is as follows [16].

For each loop

Sample a column index $j\in[k]$ uniformly. 2. 2.

Sample a row index $l\in[n]$ from distribution $\mathcal{P}_{\bm{\tilde{V}}^{\top}(:;j)}(l)=\frac{|\bm{\tilde{V}}^{\top}(l;j)|^{2}}{\|\bm{\tilde{V}}^{\top}(:;j)\|^{2}}$; 3. 3.

Compute $|\bm{\tilde{V}}(l;:)\bm{\hat{q}}_{t}|^{2}$ 4. 4.

Output $l$ with probability $\frac{|\bm{\tilde{V}}(l;:)\bm{\hat{q}}_{t}|^{2}}{\|\bm{\tilde{V}}(l;:)\|^{2}\|\bm{\hat{q}}_{t}\|^{2}}$ or restart to sample $(j,l)$ again (Back to step 1).

We execute the above loop until a sample $l$ is successfully output.

The complexity of the thin matrix-vector multiplication subroutine, namely, the required number of loops, is:

Lemma 9** ([16, Lemma 12] and [34, Proposition 6.4]).**

Let $\tilde{\bm{V}}\in\mathbb{R}^{n\times k}$ and $\bm{\hat{q}}_{t}\in\mathbb{R}^{k}$. Given $\ell_{2}$ norm sampling access to $\tilde{\bm{V}}$, we can length-square sample from $\mathcal{P}_{\tilde{\bm{V}}\bm{\hat{q}}_{t}}$ in the expected runtime complexity

[TABLE]

Appendix C General Heuristic Post-selection (Proof of Theorem 5)

Recall that the anchor $\mathcal{A}_{t}$ is defined as

[TABLE]

The goal of the general heuristic post-selection is approximating ${\mathcal{A}}_{t}$ by sampling distributions $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ with $\mathcal{O}(polylog(\max\{n_{\scalebox{0.49}{X}},n_{\scalebox{0.49}{Y}}\}))$ times, since the acquisition of the explicit form of $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ requires $\mathcal{O}(poly(n_{\scalebox{0.49}{X}},n_{\scalebox{0.49}{Y}}))$ computation complexity and collapses the desired speedup. Let $\{\bm{x}_{i}\}_{i=1}^{N_{\scalebox{0.49}{X}}}$ be $N_{\scalebox{0.49}{X}}$ examples independently sampled from $\mathcal{P}_{\bm{\hat{X}}_{t}}$ with $\mathcal{P}_{\bm{\hat{X}}_{t}}(\bm{x}=z)=|\bm{{\hat{X}}}_{t}(z)|^{2}/\|\bm{{\hat{X}}}_{t}\|^{2}$ and $N_{{\scalebox{0.49}{X}},z}$ be the number of examples taking value of $z\in[n_{\scalebox{0.49}{X}}]$ with $\sum_{z=1}^{n_{\scalebox{0.49}{X}}}N_{{\scalebox{0.49}{X}},z}=N_{{\scalebox{0.49}{X}}}$. Similarly, let $\{\bm{y}_{i}\}_{i=1}^{N_{\scalebox{0.49}{Y}}}$ be $N_{\scalebox{0.49}{Y}}$ examples independently sampled from $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ with $\mathcal{P}_{\bm{\hat{Y}}_{t}}(\bm{y}=z)=|\bm{{\hat{Y}}}_{t}(z)|^{2}/\|\bm{{\hat{Y}}}_{t}\|^{2}$ and $N_{{\scalebox{0.49}{Y}},z}$ be the number of examples taking value of $z\in[n_{\scalebox{0.49}{Y}}]$ with $\sum_{z=1}^{n_{\scalebox{0.49}{Y}}}N_{{\scalebox{0.49}{Y}},z}=N_{{\scalebox{0.49}{Y}}}$. Denote that $w_{\scalebox{0.49}{X}}$ is the total number of different indexes after sampling $\mathcal{P}_{\bm{\hat{X}}_{t}}$ with $N_{\scalebox{0.49}{X}}$ times, $I_{\scalebox{0.49}{X},v}$ and $J_{\scalebox{0.49}{X},v}$ are two indexes corresponding to the $v$-th largest value among $\{N_{\scalebox{0.49}{X},z}\}_{z=1}^{w_{\scalebox{0.49}{X}}}$ and $\mathcal{P}_{\bm{\hat{X}}_{t}}$ with $w_{\scalebox{0.49}{X}}\leq n_{\scalebox{0.49}{X}}$ and $v\in[w_{\scalebox{0.49}{X}}]$, respectively222Due to the limited sampling times, the sampled results $\{N_{\scalebox{0.49}{X},1},N_{\scalebox{0.49}{X},2},...,N_{\scalebox{0.49}{X},w_{\scalebox{0.49}{X}}}\}$ (or $\{N_{\scalebox{0.49}{Y},1},N_{\scalebox{0.49}{Y},2},...,N_{\scalebox{0.49}{Y},w_{\scalebox{0.49}{X}}}\}$) may occupy a small portion of the all $n_{\scalebox{0.49}{X}}$ (or $n_{\scalebox{0.49}{Y}}$) possible results.. Similarly, denote that $w_{\scalebox{0.49}{Y}}$ is the total number of distinguished indexes after sampling $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ with $N_{\scalebox{0.49}{Y}}$ times, the indexes $I_{\scalebox{0.49}{Y},v}$ and $J_{\scalebox{0.49}{Y},v}$ are $v$-th largest value among $\{N_{\scalebox{0.49}{Y},z}\}_{z=1}^{w_{\scalebox{0.49}{Y}}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ with $w_{\scalebox{0.49}{Y}}\leq n_{\scalebox{0.49}{Y}}$ and $v\in[w_{\scalebox{0.49}{Y}}]$, respectively. In particular, we have

[TABLE]

[TABLE]

The procedure of the general heuristic post-selection is as follows.

Sample $\mathcal{P}_{\bm{\hat{X}}_{t}}$ with $N_{\scalebox{0.49}{X}}$ times and order the sampled items to obtain $\mathcal{X}=\{I_{\scalebox{0.49}{X},1},I_{\scalebox{0.49}{X},2},...,I_{\scalebox{0.49}{X},w_{\scalebox{0.49}{X}}}\}$ with $w_{\scalebox{0.49}{X}}\leq n_{\scalebox{0.49}{X}}$; 2. 2.

Query the distribution $\mathcal{P}_{\bm{\hat{X}}_{t}}$ to obtain the value $C_{\bm{\hat{X}}_{t}}^{*}$ with $C_{\bm{\hat{X}}_{t}}^{*}=\mathcal{P}_{\bm{\hat{X}}_{t}}(\bm{x}=I_{\scalebox{0.49}{X},1})$; 3. 3.

Sample $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ with $N_{\scalebox{0.49}{Y}}$ times and order the sampled items to obtain $\mathcal{Y}=\{I_{\scalebox{0.49}{Y},1},I_{\scalebox{0.49}{Y},2},...,I_{\scalebox{0.49}{Y},w_{\scalebox{0.49}{Y}}}\}$ with $w_{\scalebox{0.49}{Y}}\leq n_{\scalebox{0.49}{Y}}$. 4. 4.

Locate $\hat{\mathcal{A}}_{t}$ with $\hat{\mathcal{A}}_{t}=I_{\scalebox{0.49}{Y},v^{*}}$ and $I_{\scalebox{0.49}{Y},v^{*}}=\arg\min_{z\in\mathcal{Y}}\left(\frac{N_{\scalebox{0.49}{Y},z}}{N_{\scalebox{0.49}{Y}}}-\hat{\xi}_{t}C_{\bm{\hat{X}}_{t}}^{*}\right)_{+}$.

An immediate observation is that, with $N_{\scalebox{0.49}{X}},~{}N_{\scalebox{0.49}{Y}}\rightarrow\infty$, we have $\hat{\mathcal{A}}_{t}=\mathcal{A}_{t}$, where

[TABLE]

The general heuristic post-selection is guaranteed by the following theorem (the formal description of Theorem 5).

Theorem 10** ((Formal) General heuristic post-selection).**

Assume that $\mathcal{P}_{\bm{\hat{X}}_{t}}$ and $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ are multinomial distributions. If $\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{\scalebox{0.49}{X},1})\geq\varepsilon_{T}$, $\left|\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{\scalebox{0.49}{X},1})-\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{\scalebox{0.49}{X},2})\right|>\varepsilon$ and $\left|\mathcal{P}_{\bm{\hat{Y}}_{t}}(J_{\scalebox{0.49}{Y},v^{*}\pm 1})-\mathcal{P}_{\bm{\hat{Y}}_{t}}(J_{\scalebox{0.49}{Y},v^{*}})\right|>\varepsilon$ for the constants $\varepsilon_{T}$ and $\varepsilon$ with $\mathcal{A}_{t}=J_{\scalebox{0.49}{Y},v^{*}}$, then for any $\delta>0$, $N_{\scalebox{0.49}{X}}\sim\mathcal{O}(\max\{1/\varepsilon,1/\varepsilon_{T}\}))$ and $N_{\scalebox{0.49}{Y}}\sim\mathcal{O}(\max\{1/\varepsilon,1/(\hat{\xi}_{t}\varepsilon_{T})\}))$, we have $\hat{\mathcal{A}}_{t}=\mathcal{A}_{t}$ with a probability at least $1-\delta$.

Remark. The physical meaning of $\varepsilon$ can be treated as the threshold of the ‘near-anchor’, that is, when the distance of a data point and the anchor point after projection is within the threshold $\varepsilon$, we say the data point can be treated as anchors. The real anchor set $\mathcal{A}$ therefore should be expanded and include these near anchors. In other words, $\varepsilon$ and $\varepsilon_{T}$ can be manually controllable. The parameter $\hat{\xi}_{t}$ is bounded as follows.

[TABLE]

In this work, we set $1/\varepsilon\sim\mathcal{O}(polylog(n_{\scalebox{0.49}{X}},n_{\scalebox{0.49}{Y}})$ and $1/\varepsilon_{T}\sim\mathcal{O}(polylog(n_{\scalebox{0.49}{X}},n_{\scalebox{0.49}{Y}})$, which gives $N_{\scalebox{0.49}{X}},~{}N_{\scalebox{0.49}{Y}}\rightarrow\mathcal{O}(\kappa_{\scalebox{0.49}{X}},polylog(\max\{n_{\scalebox{0.49}{X}},n_{\scalebox{0.49}{Y}}\}))$.

An equivalent statement of Theorem 10 is:

Problem 11**.**

How many samples, $N_{\scalebox{0.49}{X}}$ and $N_{\scalebox{0.49}{Y}}$, are required to guarantee $I_{\scalebox{0.49}{X},1}=J_{\scalebox{0.49}{X},1}$ and $I_{\scalebox{0.49}{Y},v^{*}}=J_{\scalebox{0.49}{Y},v^{*}}$, where $J_{\scalebox{0.49}{X},1}=\arg\max_{J_{{\scalebox{0.49}{X}},v}\in[n_{\scalebox{0.49}{X}}]}(\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{{\scalebox{0.49}{X}},v^{*}}))$ and $\mathcal{A}_{t}=J_{\scalebox{0.49}{Y},v^{*}}$.

We use Breteganolle-Huber-Carol inequality [36] to prove Theorem 10, i.e.,

Lemma 12** (Breteganolle-Huber-Carol inequality [36]).**

Let $\mathcal{P}_{D}$ be a multinomial distribution with $l$ event probabilities $\{\mathcal{P}_{D}(i)\}_{i=1}^{l}$. We randomly sample $N$ events from $\mathcal{P}_{D}$ and let $N_{i}$ be the number of event $i$ appeared. Then, the following holds with a probability at least $1-\delta$ for any $\delta>0$,

[TABLE]

Proof of Theorem 10.

The proof is composed of two parts. The first part is to prove that the index $I_{\scalebox{0.49}{X},1}$ with $I_{\scalebox{0.49}{X},1}=J_{\scalebox{0.49}{X},1}$ can be determined with sampling complexity $\mathcal{O}(1/(\varepsilon^{2}\varepsilon_{T}^{2}))$. The second part is to prove that the index $I_{\scalebox{0.49}{Y},w_{\scalebox{0.49}{Y}}}$ with $I_{\scalebox{0.49}{Y},w_{\scalebox{0.49}{Y}}}=\mathcal{A}_{t}$ can be determined with sampling complexity $\mathcal{O}(1/(\varepsilon^{2}\varepsilon_{T}^{2}))$.

For the first part, we split the set $\mathcal{X}$ into two subsets $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$, i.e., $\mathcal{X}_{1}=\{I_{\scalebox{0.49}{X},1}\}$ and $\mathcal{X}_{2}=\{I_{\scalebox{0.49}{X},2},...,I_{\scalebox{0.49}{X},w_{\scalebox{0.49}{X}}}\}$. The decomposition of $\mathcal{X}$ into two subsets is equivalent to setting $l=2$ in Eqn. (10). In particular, we have $N_{1}=N_{\scalebox{0.49}{X},I_{\scalebox{0.49}{X},1}}$ and $N_{2}=\sum_{v=2}^{w_{\scalebox{0.49}{X}}}N_{\scalebox{0.49}{X},I_{\scalebox{0.49}{X},v}}$. The Breteganolle-Huber-Carol inequality yields

[TABLE]

The above inequality implies that, when $\delta=4\exp{\left(\frac{-N_{\scalebox{0.49}{X}}\lambda^{2}}{2}\right)}$, we have

[TABLE]

with probability at least $1-\delta$.

The assumption $\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{\scalebox{0.49}{X},1})\geq\varepsilon_{T}$ guarantees that $J_{\scalebox{0.49}{X},1}\in\mathcal{X}$ by sampling $\mathcal{P}_{\bm{\hat{X}}_{t}}$ with $\mathcal{O}(1/\varepsilon_{T})$ times. In addition, since we have assumed that $\left|\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{\scalebox{0.49}{X},1})-\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{\scalebox{0.49}{X},2})\right|>\varepsilon$, it can be easily inferred that, when $N_{\scalebox{0.49}{X}}\geq\sqrt{\frac{2\log{(4/\delta)}}{\varepsilon^{2}}}$ with a probability at least $1-\delta$, there is one and only one value $\frac{N_{\scalebox{0.49}{X},I_{\scalebox{0.49}{X},1}}}{N_{\scalebox{0.49}{X}}}$ that is in the $\varepsilon$-neighborhood of $\mathcal{P}_{\bm{\hat{X}}_{t}}(I_{\scalebox{0.49}{X},1})$. We therefore conclude that $I_{\scalebox{0.49}{X},1}=J_{\scalebox{0.49}{X},1}$ can be guaranteed by sampling $\mathcal{P}_{\bm{\hat{X}}_{t}}$ with $N_{\scalebox{0.49}{X}}\geq\max\{\sqrt{\frac{2\log{(4/\delta)}}{\varepsilon^{2}}},\frac{1}{\varepsilon_{T}}\}$.

For the second part, we split the set $\mathcal{Y}$ into three subsets $\mathcal{Y}_{1}$, $\mathcal{Y}_{2}$ and $\mathcal{Y}_{3}$, i.e., $\mathcal{Y}_{1}=\{I_{\scalebox{0.49}{Y},1},...,I_{\scalebox{0.49}{Y},v^{*}-1}\}$, $\mathcal{Y}_{2}=\{I_{\scalebox{0.49}{Y},v^{*}}\}$, and $\mathcal{Y}_{3}=\{I_{\scalebox{0.49}{Y},v^{*}+1},...,I_{\scalebox{0.49}{Y},w_{\scalebox{0.49}{Y}}}\}$. Analogous to the above case, the decomposition of $\mathcal{Y}$ into three subsets for the case of sampling $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ is equivalent to setting $l=3$ in Eqn. (10). In particular, we have $N_{1}=\sum_{v=1}^{v^{*}}N_{\scalebox{0.49}{Y},I_{\scalebox{0.49}{Y},v}}$, $N_{2}=N_{\scalebox{0.49}{Y},I_{\scalebox{0.49}{Y},v^{*}}}$, and $N_{3}=\sum_{v=v^{*}+1}^{w_{\scalebox{0.49}{Y}}}N_{\scalebox{0.49}{Y},I_{\scalebox{0.49}{Y},v}}$. The Breteganolle-Huber-Carol inequality yields

[TABLE]

The above inequality implies that, when $\delta=4\exp{\left(\frac{-N_{\scalebox{0.49}{Y}}\lambda^{2}}{2}\right)}$, we have

[TABLE]

with probability at least $1-\delta$.

Since we have assumed $\mathcal{P}_{\bm{\hat{X}}_{t}}(J_{\scalebox{0.49}{X},1})\geq\varepsilon_{T}$, we have $\mathcal{P}_{\bm{\hat{Y}}_{t}}(J_{\scalebox{0.49}{Y},v^{*}})\geq\hat{\xi}_{t}\varepsilon_{T}$. In other words, $\mathcal{P}_{\bm{\hat{Y}}_{t}}(J_{\scalebox{0.49}{Y},v^{*}})\geq\hat{\xi}_{t}\varepsilon_{T}$ guarantees that $I_{\scalebox{0.49}{Y},v^{*}}=J_{\scalebox{0.49}{Y},v^{*}}$ by sampling $\mathcal{P}_{\bm{\hat{Y}}_{t}}$ with $\mathcal{O}(1/(\hat{\xi}_{t}\varepsilon_{T}))$ times. In addition, the assumption $\left|\mathcal{P}_{\bm{\hat{Y}}_{t}}(J_{\scalebox{0.49}{Y},v^{*}\pm 1})-\mathcal{P}_{\bm{\hat{Y}}_{t}}(J_{\scalebox{0.49}{Y},v^{*}})\right|>\varepsilon$ leads to that, when $N_{\scalebox{0.49}{Y}}\geq\max\{\sqrt{\frac{2\log{(8/\delta)}}{\varepsilon^{2}}},\frac{1}{\hat{\xi}_{t}\varepsilon_{T}}\}$ with a probability at least $1-\delta$, there is one and only one value $\frac{N_{\scalebox{0.49}{Y},I_{\scalebox{0.49}{Y},v^{*}}}}{N_{\scalebox{0.49}{Y}}}$ that is in the $\varepsilon$-neighborhood of $\mathcal{P}_{\bm{\hat{Y}}_{t}}(I_{\scalebox{0.49}{Y},v^{*}})$. We therefore conclude that $I_{\scalebox{0.49}{Y},v^{*}}=J_{\scalebox{0.49}{Y},v^{*}}$ with $J_{\scalebox{0.49}{Y},v^{*}}=\hat{\mathcal{A}}_{t}$.

Combing the results of two parts together, it can be easily inferred that, with the sampling complexity $\mathcal{\tilde{O}}(\max\{\frac{1}{\varepsilon},\frac{1}{\hat{\xi}_{t}\varepsilon_{T}}\})$, we have $\mathcal{A}_{t}=\hat{\mathcal{A}}_{t}$ with probability $1-\delta$. ∎

Appendix D Proof of Theorem 7

In this section, we give the proof of Theorem 7. We left the detailed proof of Theorem 14 and 15 in subsection D.1 and D.2, respectively.

Proof of Theorem 7.

Recall that $\bm{\hat{H}}_{t}=\tilde{\bm{V}}\bm{\hat{q}}_{\scalebox{0.49}{H,t}}$, and $\bm{\hat{q}}_{\scalebox{0.49}{H,t}}$ is an approximation of $\bm{\tilde{q}}_{\scalebox{0.49}{H,t}}=\tilde{\bm{V}}^{\top}\bm{H}_{t}$. The triangle inequality yields

[TABLE]

where $\|\mathcal{Q}\|_{TV}$ is the total variation distance of $\mathcal{Q}$. In the following, we bound the two terms on the right-hand side of Eqn. (16) respectively.

Correctness of $\|\mathcal{P}_{\bm{\tilde{H}}_{t}}-\mathcal{P}_{\bm{{H}}_{t}}\|_{TV}$**. The goal here is to prove that

[TABLE]

By Lemma 13 below, Eqn (17) follows if the following inequality holds:

[TABLE]

Finally, the inequality in Eqn. (18) is guaranteed to hold because of Theorem 14.

Lemma 13** (Lemma 6.1, [34]).**

For $\bm{x},\bm{y}\in\mathbb{R}^{n}$ satisfying $\|\bm{x}-\bm{y}\|\leq\epsilon$, the corresponding distributions $\mathcal{P}_{\bm{x}}$ and $\mathcal{P}_{\bm{y}}$ satisfy $\|\mathcal{P}_{\bm{x}},\mathcal{P}_{\bm{y}}\|_{TV}\leq\frac{2\epsilon}{\|\bm{x}\|}$.

Theorem 14**.**

Let the rank and the condition number of $\bm{H}\in\mathbb{R}^{n\times m}$ be $k$ and $\kappa$, respectively. Fix

[TABLE]

Then, Algorithm 2 yields $\|\bm{\tilde{H}}_{t}-\bm{H}_{t}\|\leq\frac{\epsilon\|\bm{H}_{t}\|}{4}$ with probability at least $(1-\eta)$.

Correctness of $\|\mathcal{P}_{\bm{\hat{H}}_{t}}-\mathcal{P}_{\bm{\tilde{H}}_{t}}\|_{TV}$**. Analogous to the above part, we bound

[TABLE]

to yield

[TABLE]

And Eqn. (19) can be obtained by the following theorem.

Theorem 15**.**

Let the rank of $\bm{H}\in\mathbb{R}^{n\times m}$ be $k$. Set the number of samplings in the inner product subroutine as

[TABLE]

Then, Algorithm 2 yields $\|\bm{\hat{H}}_{t}-\bm{\tilde{H}}_{t}\|\leq\frac{\epsilon}{4}\|\bm{\tilde{H}}_{t}\|$ with at least $1-\delta$ success probability.

Finally, Eqn. (16) holds combining Eqn. (17) with Eqn. (20). ∎

D.1 Proof of Theorem 14

Due to $\bm{H}_{t}=\bm{V}\bm{V}^{\top}\bm{H}\bm{B}_{t}$, we have

[TABLE]

where $\Pi_{(\bm{H})}=\sum_{i}\bm{v}^{(i)}\bm{v}^{(i)\top}$ and $\bm{v}^{(i)}$ is the left singular vectors of $\bm{H}$. The first inequality of Eqn. (21) is obtained by exploiting the submultiplicative property of spectral norm [18], i.e., for any matrix $\bm{M}\in\mathbb{R}^{n\times m}$ and any vector $\bm{z}\in\mathbb{R}^{m}$, we have $\|\bm{M}\bm{z}\|\leq\|\bm{M}\|_{2}\|\bm{z}\|$. The second inequality of Eqn. (21) comes from the submultiplicative property of spectral norm and $\|\bm{B}_{t}\|=1$. To achieve $\|\bm{\tilde{H}}_{t}-\bm{H}_{t}\|\leq\frac{\epsilon}{4}\|\bm{H}_{t}\|$ in Eqn. (18), Eqn. (21) indicates that the approximated left singular matrix $\tilde{\bm{V}}$ should satisfy

[TABLE]

The spectral norm $\|\sum_{i=1}^{k}\tilde{\bm{v}}^{(i)}\tilde{\bm{v}}^{(i)\top}-\Pi_{(\bm{H})}\|_{2}$ can be quantified as:

Theorem 16**.**

Suppose that the rank of $\bm{H}$ is $k$, and $\tilde{\bm{v}}^{(i)}$ refers to a approximated singular vector of $\bm{H}$ such that

[TABLE]

Then, we have

[TABLE]

The proof of Theorem 16 is given in Subsection D.1.1.

Theorem 16 implies that to achieve Eqn. (22) (or equivalently, Eqn. (18)), we should bound $\alpha$ as

[TABLE]

We use the following lemma to give an explicit representation of $\alpha$ by the sampled matrix $\bm{R}$ and $\bm{C}$,

Lemma 17**.**

Suppose that $\bm{{\omega}}^{(l)}$ refers to the right singular vector of $\bm{C}$ such that $\Pi_{(\bm{C})}=\sum_{l}\bm{\omega}^{(l)}\bm{\omega}^{(l)\top}$ and $\bm{\bm{\omega}}^{(i)\top}\bm{C}^{\top}\bm{C}\bm{\omega}^{(j)}=\delta_{ij}(\sigma^{(i)})^{2}$, where $(\sigma^{(i)})^{2}\geq 4/(5\kappa^{2})$. Suppose that the rank of both $R$ and $C$ is $k$ and

[TABLE]

Let $\tilde{\bm{v}}^{(l)}:=\bm{R}\bm{\omega}^{(l)}/\sigma^{(l)}$, then we have

[TABLE]

The proof of Lemma 17 is presented in Subsection D.1.2.

In conjunction with Eqn. (23) and Eqn. (26), we set $\alpha={(5\gamma\kappa^{2})}/{4}$ and rewrite Eqn. (25) as

[TABLE]

In other words, when $\gamma\leq\frac{3\epsilon}{85k\kappa^{2}}$, Eqn. (18) is achieved so that $\|\bm{\tilde{H}}_{t}-\bm{H}_{t}\|\leq\frac{\epsilon}{4}\|\bm{H}_{t}\|$. Recall that $\gamma$ is quantified by $\|\bm{R}^{\top}\bm{R}-\bm{C}\bm{C}^{\top}\|_{2}$ as defined in Eqn. (26), we use the following theorem to bound $\gamma$, i.e.,

Theorem 18**.**

Given a nonnegative matrix $\bm{H}\in\mathbb{R}^{n\times m}$, let $\bm{R}\in\mathbb{R}^{n\times s}$, $\bm{C}\in\mathbb{R}^{s\times s}$ be the sampled matrix following Algorithm 1. Setting $s$ as $s=\frac{85^{2}k^{2}\kappa^{4}\ln(8n/\eta)\|\bm{H}\|_{F}^{2}}{9\epsilon^{2}}$, with probability at least $(1-\eta)$, we always have $\|\bm{R}^{\top}\bm{R}-\bm{C}^{\top}\bm{C}\|\leq\gamma$,

The proof of Theorem 18 is given in D.1.3.

Combining the result of Theorem 16 and Lemma 17, we know that with sampling $s$ rows of $\bm{H}$, the approximated distribution is $\frac{\epsilon}{2}$-close to the desired result, i.e.,

[TABLE]

D.1.1 Proof of Theorem 16

We first introduce a lemma to facilitate the proof of Theorem 16, i.e.,

Lemma 19** (Adapted from Lemma 5, [16]).**

Let $\bm{A}$ be a matrix of rank at most $k$, and suppose that $\bm{W}$ has $k$ columns that span the row and column spaces of $\bm{A}$. Then $\|\bm{A}\|_{2}\leq\|(\bm{W}^{\top}\bm{W})^{-1}\|_{2}\|\bm{W}^{\top}\bm{AW}\|_{2}$ .

Proof of Theorem 16.

The main procedure to prove this theorem is as follows. By employing the Lemma 19, we can set $\bm{A}$ and $\bm{W}$ as

[TABLE]

and then bound $\|\bm{W}^{\top}\bm{AW}\|_{2}$ and $\|(\bm{W}^{\top}\bm{W})^{-1}\|_{2}$ separately. Lastly, we combine the two results to obtain the bound $\|\sum_{i=1}^{k}\tilde{\bm{v}}^{(i)}\tilde{\bm{v}}^{(i)\top}-\Pi_{(\bm{H})}\|_{2}\leq\frac{17k\alpha}{3}$ in Eqn. (24).

Following the above observation, we first bound the term $\|\bm{W}^{\top}\bm{AW}\|_{2}$. We rewrite $\bm{W}^{\top}\bm{AW}$ as

[TABLE]

The entry $\bm{A}(i,j)$ of $\bm{A}$ with $\bm{A}(i,j)=(\tilde{\bm{v}}^{(i)\top}A\tilde{\bm{v}}^{(j)})$ is bounded by $\alpha$, i.e.,

[TABLE]

The first equivalence of Eqn. (D.1.1) comes from the definition of $\bm{A}$, and the second equivalence employs $\tilde{\bm{v}}^{(i)\top}\Pi_{(\bm{H})}\tilde{\bm{v}}^{(j)}=\tilde{\bm{v}}^{(i)\top}\tilde{\bm{v}}^{(j)}.$ The first inequality of Eqn. (D.1.1) exploits triangle inequality and Eqn. (23), i.e.,

[TABLE]

The second inequality of Eqn. (D.1.1) directly comes from the triangle inequality. The last second inequality of Eqn. (D.1.1) employs the inequality $\left|\sum_{t^{\prime}=\{i,j\}}\tilde{\bm{v}}^{(i)\top}\tilde{\bm{v}}^{(t^{\prime})}\tilde{\bm{v}}^{(t^{\prime})\top}\tilde{\bm{v}}^{(j)}-\delta_{ij}\right|\leq 3\alpha$ for both the case $i=j$ and $i\neq j$, guaranteed by Eqn. (23) and $\alpha^{2}<\alpha$. Specifically, for the case $i\neq j$, we bound $\left|\sum_{t^{\prime}=\{i,j\}}\tilde{\bm{v}}^{(i)\top}\tilde{\bm{v}}^{(t^{\prime})}\tilde{\bm{v}}^{(t^{\prime})\top}\tilde{\bm{v}}^{(j)}-\delta_{ij}\right|$ as

[TABLE]

For the case $i=j$, we bound $\left|\sum_{t^{\prime}=\{i,j\}}\tilde{\bm{v}}^{(i)\top}\tilde{\bm{v}}^{(t^{\prime})}\tilde{\bm{v}}^{(t^{\prime})\top}\tilde{\bm{v}}^{(j)}-\delta_{ij}\right|$ as

[TABLE]

The last inequality of Eqn. (D.1.1) comes from

[TABLE]

since

[TABLE]

Combing Eqn. (28) and Eqn. (D.1.1), we immediately have

[TABLE]

In addition, from Eqn. (23) and the definition of $\bm{W}$, we can obtain

[TABLE]

and then $\left\|(\bm{W}^{\top}\bm{W})^{-1}\right\|_{2}\leq 4/3$. Combining the two result, we have

[TABLE]

∎

D.1.2 Proof of Lemma 17

Proof of Lemma 17.

The inequality $\left|\tilde{\bm{v}}^{(i)\top}\tilde{\bm{v}}^{(j)}-\delta_{ij}\right|\leq\frac{5\kappa^{2}\gamma}{4}$ in Eqn. (26) can be proved following the definition of $\tilde{\bm{v}}^{(i)}$. Mathematically, we have

[TABLE]

The first equivalence of Eqn. (D.1.2) comes from the definition of $\tilde{\bm{v}}^{(i)}$. The first inequality of Eqn. (D.1.2) is derived by employing $\|\bm{R}^{\top}\bm{R}-\bm{C}^{\top}\bm{C}\|\leq\gamma$, i.e.,

[TABLE]

The last second equivalence of Eqn. (D.1.2) employs

[TABLE]

The last inequality of Eqn. (D.1.2) uses

[TABLE]

∎

D.1.3 Proof of Theorem 18

We introduce the following lemma to facilitate the proof.

Lemma 20** (Adapted from Theorem 4.4, [19]).**

Given any matrix $\bm{R}\in\mathbb{R}^{n\times s}$. Let ${\bm{C}}\in\mathbb{R}^{s\times s}$ be obtained by length-squared sampling with $\mathbb{E}(\bm{C}^{\top}\bm{C})={\bm{R}}^{\top}{\bm{R}}$. Then, for all $\epsilon\in[0,\|\bm{R}\|_{2}/\|\bm{R}\|_{F}]$, we have

[TABLE]

Hence, for $s\geq 4\ln{(2n/\eta)/\epsilon^{2}}$, with probability at least $(1-\eta)$ we have $\left\|\bm{C}^{\top}\bm{C}-{\bm{R}}^{\top}{\bm{R}}\right\|_{2}\leq\epsilon\|\bm{R}\|\|\bm{R}\|_{F}$.

Proof of Theorem 18.

The Lemma 20 indicates that the sample complexity of $s$ determines $\|\bm{C}^{\top}\bm{C}-{\bm{R}}^{\top}{\bm{R}}\|_{2}$. With setting $\gamma=\epsilon\|\bm{R}\|_{2}\|\bm{R}\|_{F}$, we have

[TABLE]

Let the right hand side of the above inequality be $\eta$, i.e.,

[TABLE]

where the inequality comes from $\|\bm{R}\|_{F}\leq\|\bm{R}\|_{F}$ and $\|\bm{R}\|_{F}=\|\bm{H}\|_{F}$. Therefore, with setting $s$ as

[TABLE]

we have $\|{\bm{R}}^{\top}{\bm{R}}-\bm{C}^{\top}\bm{C}\|_{2}\leq\gamma$ with probability at least $(1-\eta)$. ∎

D.2 Proof of Theorem 15

Proof of Theorem 15.

We first give the upper bound of the term $\left\|\hat{\bm{H}}_{t}-\tilde{\bm{H}}_{t}\right\|$, i.e.,

[TABLE]

The first inequality comes from the the submultiplicative property of spectral norm. The second inequality supported by Eqn. (27) with

[TABLE]

Following the definition of $\ell_{2}$ norm, we have

[TABLE]

Denote the additive error

[TABLE]

we rewrite Eqn. (37) as

[TABLE]

An observation of the above equation is that we have $\left\|\hat{\bm{H}}_{t}-\tilde{\bm{H}}_{t}\right\|\leq\epsilon/4$ if

[TABLE]

We use the result of the inner product subroutine to quantify the required number of samplings to achieve Eqn. (39). The conclusion of Lemma 8 is that when

[TABLE]

we have $|\bm{\tilde{q}}_{{\scalebox{0.49}{H,t}}}(i)-\bm{\hat{q}}_{{\scalebox{0.49}{H,t}}}(i)|\leq\epsilon$ with at least $1-\delta$ success probability. With substituting $\epsilon$ by $\epsilon^{\prime}$, we immediately obtain $|\bm{\tilde{q}}_{{\scalebox{0.49}{H,t}}}(i)-\bm{\hat{q}}_{{\scalebox{0.49}{H,t}}}(i)|\leq\epsilon^{\prime}$, where the required number of samplings is

[TABLE]

∎

Appendix E The Complexity of The Algorithm (Proof of Theorem 6)

Proof of Theorem 6.

As analyzed in the main text, the complexity of the proposed algorithm is dominated by four operations in the preprocessing step and the divide step, i.e., finding the left singular vectors $\tilde{\bm{V}}$, estimating the inner product to build $\hat{\bm{q}}_{t}$, preparing the approximated probability distribution $\mathcal{P}_{\hat{\bm{H}}}$, and estimating the rescale factor $\hat{\xi}_{t}$. We evaluate the computation complexity of these four operations separately and then give the overall computation complexity of our algorithm.

In this subsection, we first evaluate the computation complexity of these four parts separately and then combine the results to give the computation complexity of our algorithm. Due to same reconstruction rule, we use a general setting $\bm{H}\in\mathbb{R}^{n\times m}$ that can either be $\bm{X}$ or $\bm{Y}$ to evaluate the computation complexity for the four parts.

Complexity of Finding $\tilde{\bm{V}}$. Supported by the $\ell_{2}$ norm sampling operations, the matrix $\bm{C}$ can be efficiently constructed following Algorithm 1, where $\mathcal{O}(2s\log^{2}(mn))$ query complexity is sufficient. Applying SVD onto $\bm{C}\in\mathbb{R}^{s\times s}$ with $s=\frac{85^{2}k^{2}\kappa^{4}\ln(8n/\eta)\|\bm{H}\|_{F}^{2}}{9\epsilon^{2}}$ generally costs

[TABLE]

runtime complexity. Once we obtain such the SVD result of $\bm{C}$, the approximated left singular vectors $\tilde{\bm{V}}$ can be implicitly represented, guaranteed by the following Lemma:

Lemma 21** (Adapted from [34]).**

Let the given dataset support the $\ell_{2}$ norm sampling operations along with the description of $\tilde{\bm{V}}\in\mathbb{R}^{n\times s}$, We can sample from any $\tilde{\bm{v}}^{(t)}$ in $\mathcal{O}(Ks^{2})$ expected queries with $K=\kappa\|\bm{H}\|_{F}^{2}$ and query for any particular entry $\tilde{\bm{V}}(i,j)$ in $\mathcal{O}(s)$ queries.

Complexity of Estimating $\bm{\tilde{q}}_{t}$ by $\bm{\hat{q}}_{t}$. The runtime complexity to estimate $\bm{\tilde{q}}_{t}$ by $\bm{\hat{q}}_{t}$ obeys the following corollary, i.e.,

Corollary 22**.**

Let $\bm{B}_{t}\in\mathbb{R}^{m\times 1}$ be the input vector, $\bm{H}\in\mathbb{R}^{n\times m}$ be the input matrix with rank $k$, and $\bm{\tilde{V}}\in\mathbb{R}^{n\times k}$ be the approximated left singular matrix. We can estimate $\tilde{\bm{q}}_{t}=\tilde{\bm{V}}^{\top}\bm{H}\bm{B}_{t}$ by $\hat{\bm{q}}_{t}$ to precision $\epsilon$ with probability at least $1-\delta$ using

[TABLE]

runtime complexity.

Proof of Corollary 22.

The proof of Corollary 22 employs the result of Theorem 15 and Lemma 8.

Theorem 15 indicates that, to estimate $\tilde{\bm{q}}_{t}(i)$ by $\hat{\bm{q}}_{t}(i)$, the required number of samplings is

[TABLE]

Following the result of Lemma 8, the runtime complexity to obtain $\hat{\bm{q}}_{t}(i)$ is

[TABLE]

Since $\|\tilde{\bm{v}}^{(i)}\|\leq\sqrt{1+\epsilon}$ for any $\tilde{\bm{v}}^{(i)}$ indicated by Eqn. (37), we rewrite Eqn. (41) as

[TABLE]

We now quantify the access cost of $\bm{H}$, $\bm{B}_{t}$ and $\bm{\tilde{V}}$ to give an explicit bound of Eqn. (41). We have $L(\bm{H})=\mathcal{O}(\log^{2}(nm))$ and $Q(\bm{B}_{t})=\mathcal{O}(\log(m))$, since $\bm{H}$ and $\bm{B}_{t}$ are stored in BNS data structure. We have $Q(\bm{\tilde{V}})=\mathcal{O}(s)$ supported by Lemma 21. Combing the above access cost and Eqn. (41), the runtime complexity to estimate $\hat{\bm{q}}_{t}(i)$ is

[TABLE]

Since each entry can be computed in parallel, the runtime complexity to obtain $\bm{\hat{q}}_{t}$ is also

[TABLE]

∎

Complexity of Sampling from $\mathcal{P}_{\bm{\hat{H}}_{t}}$. Recall that the definition of ${\bm{\tilde{H}}_{t}}$ is ${\bm{\hat{H}}_{t}}=\tilde{V}\bm{\hat{q}}_{t}$. We first evaluate the computation complexity to obtain one sample from $\mathcal{P}_{\bm{\hat{H}}_{t}}$. From Lemma 9, we know the expected runtime complexity to sample from $\mathcal{P}_{\bm{\hat{H}}_{t}}$ is $\mathcal{O}(\frac{k\|\bm{\hat{q}}_{t}\|^{2}}{\|\tilde{\bm{V}}\bm{\hat{q}}_{t}\|^{2}}(S(\tilde{\bm{V}})+kQ(\tilde{\bm{V}})))$. Specifically, we have

[TABLE]

where the first inequality employs the triangle inequality, the second inequality employs the submultiplicative property of spectral norm, and the third inequality utilizes $\|\bm{H}_{t}\|\leq\|\bm{H}\|_{2}\|\bm{B}_{t}\|\leq\|\bm{H}\|_{F}\|\bm{B}_{t}\|$ with $\|\bm{B}_{t}\|=1$. Concurrently, we have $\|\tilde{\bm{V}}\bm{q}_{t}\|=\Omega(1)$. Employing Lemma 21 to quantify $S(\tilde{\bm{V}})$ and $Q(\tilde{\bm{V}})$, the complexity to obtain a sample from $\mathcal{P}_{\bm{\hat{H}}_{t}}$ is

[TABLE]

With substituting $s$ with its explicit representation in Theorem 7, the complexity is

[TABLE]

Following the result of the general heuristic post-selection method in Theorem 10, we sample the distribution $\mathcal{P}_{\bm{\hat{H}}_{t}}$ with $N\sim\mathcal{\tilde{O}}(\frac{1}{\varepsilon})$ times in parallel, which gives the runtime complexity

[TABLE]

Complexity of estimating $\xi_{t}$ by $\hat{\xi}_{t}$. For the general case with $\bm{X}\neq\bm{Y}$, we should estimate the $\ell_{2}$ norm of their project results, i.e., $\|\bm{H}_{t}\|^{2}=\bm{H}_{t}^{\top}\bm{H}_{t}$ that $\bm{H}$ can either be $\bm{X}$ or $\bm{Y}$. Recall that the explicit representation of the approximated result is $\|\bm{\hat{H}}_{t}\|$ with $\|\bm{H}_{t}\|=\bm{\hat{q}}_{t}^{\top}\bm{\tilde{V}}^{\top}\bm{\tilde{V}}\bm{\hat{q}}_{t}$, where the $(j,j)$-th entry of $\bm{\tilde{V}}^{\top}\bm{\tilde{V}}$ is $\|\tilde{\bm{v}}^{(j)}\|^{2}$ and the else entries are zero. An immediate observation is that $\|\bm{H}_{t}\|$ can be obtained by using the inner product subroutine in Lemma 8.

We first calculate the sample and query complexity to query the $(j,j)$-th entry of $\bm{\tilde{V}}^{\top}\bm{\tilde{V}}$, namely, the query and sample complexity to obtain the inner product of $\tilde{\bm{v}}^{(j)}$. By employing the result of the inner product subroutine, with removing $\bm{H}$ and setting both $\bm{B}_{t}$ and $\bm{\tilde{v}}^{(i)}$ as $\bm{\tilde{v}}^{(j)}$, we can estimate $\tilde{\bm{v}}^{(j)\top}\tilde{\bm{v}}^{(j)}$ to precision $\epsilon$ with probability at least $1-\delta$ in time

[TABLE]

where we use Lemma 21 to get $Q(\tilde{\bm{v}}^{(i)})=\mathcal{O}(s)$ and the inequality comes from $\|\tilde{\bm{v}}^{(i)}\|\leq(1+\epsilon)$ in Eqn. (21). We store $k$ nonzero entries of $\bm{\tilde{V}}^{\top}\bm{\tilde{V}}$ in memory.

We next use the inner product subroutine to obtain $\|\bm{\hat{H}}_{t}\|$ with $\|\bm{\hat{H}}_{t}\|=\bm{\hat{q}}_{t}^{\top}\bm{\tilde{V}}^{\top}\bm{\tilde{V}}\bm{\hat{q}}_{t}$. Since both $\bm{\hat{q}}_{t}$ and $\bm{\tilde{V}}^{\top}\bm{\tilde{V}}$ are stored in memory, we have $L(\bm{\tilde{V}}^{\top}\bm{\tilde{V}})=\mathcal{O}(k)$ and $Q(\bm{\hat{q}}_{t})=\mathcal{O}(1)$. Following the result of Lemma 8, we can estimate $\|\hat{\bm{H}}_{t}\|$ to precision $\epsilon$ with probability at least $1-\delta$ with runtime complexity

[TABLE]

where the inequality comes from $\|\tilde{\bm{v}}^{(i)}\|\leq(1+\epsilon)$ in Eqn. (21) and Eqn. (43).

Combining Eqn. (44) and Eqn. (E), the computation complexity to obtain the $\ell_{2}$ norm of $\|\hat{\bm{H}}_{t}\|$ is

[TABLE]

The overall complexity of our algorithm. An immediate observation of the above four parts is that the query complexity and runtime complexity of our algorithm is denominated by the complexity of finding $\tilde{V}$, i.e., $\mathcal{O}({s^{3}})=\tilde{\mathcal{O}}(\frac{85^{6}k^{6}\kappa^{12}\|\bm{H}\|_{F}^{6}}{9^{3}\epsilon^{6}})$ . Since $\bm{H}$ is a general setting that can either be $\bm{X}$ or $\bm{Y}$, the runtime complexity for our algorithm is

[TABLE]

∎