Painless Breakups -- Efficient Demixing of Low Rank Matrices

TL;DR

This paper introduces two efficient algorithms for demixing multiple low-rank matrices from combined measurements, with theoretical guarantees and applications in quantum tomography and IoT.

Contribution

The paper proposes novel hard thresholding algorithms for low rank matrix demixing with proven linear convergence guarantees.

Findings

Algorithms successfully demix low-rank matrices in simulations.

Theoretical analysis confirms linear convergence under certain conditions.

Applications demonstrated in quantum tomography and IoT contexts.

Abstract

Assume we are given a sum of linear measurements of $s$ different rank-$r$ matrices of the form $y = \sum_{k=1}^{s} \mathcal{A}_k ({X}_k)$. When and under which conditions is it possible to extract (demix) the individual matrices ${X}_k$ from the single measurement vector ${y}$? And can we do the demixing numerically efficiently? We present two computationally efficient algorithms based on hard thresholding to solve this low rank demixing problem. We prove that under suitable conditions these algorithms are guaranteed to converge to the correct solution at a linear rate. We discuss applications in connection with quantum tomography and the Internet-of-Things. Numerical simulations demonstrate empirically the performance of the proposed algorithms.

Figures

Equations

$$\mathbf{y}=\sum _{k=1}^s{\mathcal{A}}_k({\mathbf{X}}_k),$$

$$\mathbf{y}=\mathbf{A}\left(\sum _{k=1}^s{\mathbf{U}}_k{\mathbf{X}}_k+{\mathbf{w}}_k\right),$$

$$r$$

$${\mathcal{A}}_k(\mathbf{Z})=\frac{1}{\sqrt{m}}\left[\begin{array}{c}⟨{\mathbf{A}}_{k,1},\mathbf{Z}⟩\\ ⋮\\ ⟨{\mathbf{A}}_{k,m},\mathbf{Z}⟩\end{array}\right].$$

$$(1-{\delta }_r)\sum _{k=1}^s{\Vert {\mathbf{Z}}_k\Vert }_F^2\le {\Vert \sum _{k=1}^s{\mathcal{A}}_k({\mathbf{Z}}_k)\Vert }^2\le (1+{\delta }_r)\sum _{k=1}^s{\Vert {\mathbf{Z}}_k\Vert }_F^2$$

$$(1-{\delta }_r){\Vert {\mathbf{Z}}_k\Vert }_F^2\le {\Vert {\mathcal{A}}_k({\mathbf{Z}}_k)\Vert }^2\le (1+{\delta }_r){\Vert {\mathbf{Z}}_k\Vert }_F^2$$

$$m\ge C{\delta }^{-2}(2n+1)rs\log (1/\delta ),$$

$$\underset{\{{\mathbf{Z}}_k{\}}_{k=1}^s}{\min }{\Vert \mathbf{y}-\sum _{k=1}^s{\mathcal{A}}_k({\mathbf{Z}}_k)\Vert }^2\text{mbox}subjectto\mathrm{rank}({\mathbf{Z}}_k)\le r.$$

$$T_{k,l}=\{{\mathbf{U}}_{k,l}{\mathbf{Z}}_1^{\ast }+{\mathbf{Z}}_2{\mathbf{V}}_{k,l}^{\ast }|{\mathbf{Z}}_1\in {\mathbb{R}}^{n\times r},{\mathbf{Z}}_2\in {\mathbb{R}}^{n\times r}\}.$$

$${\gamma }_1=\frac{4{\delta }_{3r}}{1-{\delta }_{3r}}.$$

$${\Vert \left[\begin{array}{c}{\mathbf{X}}_{1,l+1}-{\mathbf{X}}_1\\ ⋮\\ {\mathbf{X}}_{s,l+1}-{\mathbf{X}}_s\end{array}\right]\Vert }_F\le {\gamma }_1{\Vert \left[\begin{array}{c}{\mathbf{X}}_{1,l}-{\mathbf{X}}_1\\ ⋮\\ {\mathbf{X}}_{s,l}-{\mathbf{X}}_s\end{array}\right]\Vert }_F.$$

$${\Vert \left[\begin{array}{c}{\mathbf{X}}_{1,l}-{\mathbf{X}}_1\\ ⋮\\ {\mathbf{X}}_{s,l}-{\mathbf{X}}_s\end{array}\right]\Vert }_F\le {\gamma }_1^l{\Vert \left[\begin{array}{c}{\mathbf{X}}_{1,0}-{\mathbf{X}}_1\\ ⋮\\ {\mathbf{X}}_{s,0}-{\mathbf{X}}_s\end{array}\right]\Vert }_F+\frac{\xi }{1-{\gamma }_1}\Vert \mathbf{e}\Vert ,$$

$$\xi =\frac{2\sqrt{1+{\delta }_{3r}}}{1-{\delta }_{2r}}.$$

$${\mathbf{W}}_{k,l}=\left[\begin{array}{cc}{\mathbf{U}}_{k,l}& {\mathbf{Q}}_1\end{array}\right]{\mathbf{M}}_{k,l}\left[\begin{array}{c}{\mathbf{V}}_{k,l}^{\ast }\\ {\mathbf{Q}}_2^{\ast }\end{array}\right]$$

$${\mathbf{M}}_{k,l}=\left[\begin{array}{cc}{\mathbf{\Sigma }}_{k,l}+{\mathbf{\xi }}_{k,l}& {\mathbf{R}}_1^{\ast }\\ {\mathbf{R}}_2^{\ast }& 0\end{array}\right].$$

$${\gamma }_2=2\left(\frac{2{\delta }_{2r}}{1-{\delta }_{2r}}+\frac{{\delta }_{3r}}{1-{\delta }_{2r}}+2{\delta }_{3r}\sqrt{rs}\frac{{\sigma }_{\max }}{{\sigma }_{\min }}\right),$$

$${\Vert \left[\begin{array}{c}{\mathbf{X}}_{1,l+1}-{\mathbf{X}}_1\\ ⋮\\ {\mathbf{X}}_{s,l+1}-{\mathbf{X}}_s\end{array}\right]\Vert }_F\le {\gamma }_2{\Vert \left[\begin{array}{c}{\mathbf{X}}_{1,l}-{\mathbf{X}}_1\\ ⋮\\ {\mathbf{X}}_{s,l}-{\mathbf{X}}_s\end{array}\right]\Vert }_F.$$

$${\delta }_{3r}\lesssim \frac{{\sigma }_{\min }}{{\sigma }_{\max }}\frac{1}{\sqrt{rs}}.$$

$$\frac{\Vert \mathbf{y}-{\sum }_{k=1}^s{\mathcal{A}}_k({\mathbf{X}}_{k,l})\Vert }{\Vert \mathbf{y}\Vert }\le 10^{-4}.$$

$$\frac{\sqrt{{\sum }_{k=1}^s{\Vert {\mathbf{X}}_{k,l}-{\mathbf{X}}_k\Vert }_F^2}}{\sqrt{{\sum }_{k=1}^s{\Vert {\mathbf{X}}_k\Vert }_F^2}}\le 10^{-2}.$$

$$T_{k,l}=\{{\mathbf{U}}_{k,l}{\mathbf{Z}}^{\ast }+\mathbf{Z}{\mathbf{U}}_{k,l}^{\ast }|\mathbf{Z}\in {\mathbb{C}}^{n\times r}\}.$$

$${\mathbf{W}}_{k,l}=\left[\begin{array}{cc}{\mathbf{U}}_{k,l}& \mathbf{Q}\end{array}\right]{\mathbf{M}}_{k,l}\left[\begin{array}{c}{\mathbf{U}}_{k,l}^{\ast }\\ {\mathbf{Q}}^{\ast }\end{array}\right],$$

$${\mathbf{M}}_{k,l}=\left[\begin{array}{cc}{\mathbf{\Lambda }}_{k,l}+{\mathbf{\xi }}_{k,l}& {\mathbf{R}}^{\ast }\\ {\mathbf{R}}^{\ast }& 0\end{array}\right].$$

$$\Pi =\{\mathbf{Z}\in {\mathbb{C}}^{n\times n}|\mathbf{Z}={\mathbf{Z}}^{\ast },\mathrm{rank}(\mathbf{Z})\le r,\text{mbox}and\mathbf{Z}⪰0\};$$

$$\Delta =\left\{\mathbf{x}\in {\mathbb{R}}^r|\sum _{i=1}^r{x}_i=1,x_i\ge 0\right\}.$$

$${\mathbf{A}}_{k,p}=2^{-q/2}{\otimes }_{i=1}^q{\sigma }_a,$$

$${\sigma }_1=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],{\sigma }_2=\left[\begin{array}{cc}0& -i\\ i& 0\end{array}\right],{\sigma }_3=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],{\sigma }_4=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].$$

$$\mathbf{y}=\sum _k{\mathbf{z}}_k\ast {\mathbf{g}}_k.$$

$$\frac{1}{\sqrt{m}}\hat{\mathbf{y}}=(\mathbf{F}{\mathbf{C}}_k{\mathbf{x}}_k)\odot \mathbf{B}{\mathbf{h}}_k,$$

$$[(\mathbf{F}{\mathbf{C}}_k{\mathbf{x}}_k)\odot \mathbf{B}{\mathbf{h}}_k{]}_p=({\mathbf{c}}_{k,p}^{\ast }{\mathbf{x}}_k)\cdot ({\mathbf{b}}_p^{\ast }{\mathbf{h}}_k)=⟨{\mathbf{c}}_{k,p}{\bar{\mathbf{b}}}_p^{\ast },{\mathbf{x}}_k{\bar{\mathbf{h}}}_k^{\ast }⟩=⟨{\mathbf{c}}_{k,p}{\bar{\mathbf{b}}}_p^{\ast },{\mathbf{X}}_k⟩,$$

Full text

1 Introduction

Demixing problems appear in a wide range of areas, including audio source separation [16], image processing [3], communications engineering [25], and astronomy [20]. A vast amount of recent literature focuses on the demixing of signals with very different properties. For instance one may be interested in the demixing of a signal that is the sum of spikes and sines [17] or the separation of a matrix that is the sum of a sparse matrix and a low rank matrix [4]. This paper focuses on the demixing problem where the signals to be separated all have similar properties, i.e., they are all low rank matrices. Assume we are given a sum of linear measurements

[TABLE]

where $\{\mathcal{A}_{k}\}_{k=1}^{s}$ is a set of linear operators from $n\times n$ matrices to $m$-dimensional vectors and the constituents $\{\bm{X}_{k}\}_{k=1}^{s}$ are unknown rank-$r$ matrices. Our goal is to extract each constituent matrix $\bm{X}_{k}$ from the single observation vector $\bm{y}$. We face different challenges in this homogeneous scenario, and one way to make the separation possible is by ensuring that the sensing operators $\mathcal{A}_{k}$ are sufficiently different from each other. The problem of demixing a sum of low rank matrices arises in quantum tomography [8], dictionary learning [12] and wireless communications [29].

In quantum tomography one tries to reconstruct an unknown quantum state from experimental data [13]. The state of a quantum system in quantum mechanics is often described by a density matrix, which is a positive semi-definite Hermitian matrice with unit trace. Many density matrices of interest are (or at least approximately) low rank—for instance pure states can be represented by rank-one matrices. Hence, in quantum tomography the goal is to reconstruct low rank matrices from a set of linear measurements, and ideally by using as few measurements as possible. In [8], a specific measurement protocol for quantum tomography is described in which states get mixed together. In this case one has to reconstruct and demix quantum states (i.e., low rank matrices) from a set of linear measurements.

Demixing problems of the form (1) are also expected to arise in the future Internet-of-Things (IoT). The IoT will connect billions of wireless devices, which is far more than the current wireless systems can technically and economically accommodate. One of the many challenges in the design of the IoT will be its ability to manage the massive number of sporadic traffic generating devices which are inactive most of the time, but regularly access the network for minor updates with no human interaction [29]. It is common understanding among communication engineers that this traffic cannot be handled within the current random access procedures. Dimensioning the channel access according to classical information and communication theory results in a severe waste of resources which does not scale towards the requirements of the IoT. This means among others that the overhead caused by the exchange of certain types of information between transmitter and receiver, such as channel estimation, assignment of data slots, etc, has to be avoided as much as possible. Without explicit channel estimation the receiver needs to blindly deconvolve (to undo the effect of the channel) and demix the signals arriving from many different devices. We will describe in Section 3 how this blind deconvolution-demixing problem can be phrased as a demixing problem of rank-one matrices.

1.1 State of the Art

The low rank matrix demixing problem in (1) is a generalization of the well-known low rank matrix recovery problem [19]. The task in low rank matrix recovery is to reconstruct a single low rank matrix $\bm{X}$ from a few linear measurements $\bm{y}=\mathcal{A}(\bm{X})$; that is, we have $s=1$ in (1). This problem finds applications in a wide range of disciplines, including quantum tomography [13], and image processing [19, 30]. Many different algorithms have been proposed for its solution [30], including convex optimization based methods [6] and non-convex methods such as thresholding algorithms [2, 22]. The thresholding-based algorithms proposed in this paper can be viewed as natural generalizations of the latter algorithms.

There is a plethora of literature of various demixing problems, such as  demixing of a low rank matrix and a sparse matrix, see for example [4]. However, most of these problems differ considerably from the one studied in this paper, since they focus on the separation of two signals with complementary properties. Therefore we will not discuss them here in any further detail.

Two of the first papers that consider the demixing of a sum of low rank matrices are [28] and [18]. In [18], the authors consider a general demixing framework of the form

[TABLE]

where $\bm{A}$ is a right-invertible matrix, the matrices $\{\bm{U}_{k}\}_{k=1}^{s}$ are random unitary matrices, the matrices $\{\bm{X}_{k}\}_{k=1}^{s}$ are the signals of interest and the vectors $\{\bm{w}_{k}\}_{k=1}^{s}$ are noise. In the framework by McCoy and Tropp the signals $\bm{X}_{k}$ are assumed to be highly structured, which includes the scenario where all the $\bm{X}_{k}$ are rank-$r$ matrices. Thus, in that case the setup (2) is a special instance of (1). The focus of [18] is on theoretical performance bounds of convex approaches for solving (2). In a nutshell the authors derive compelling quantitative phase transition bounds which show that demixing can succeed if and only if the dimension of the observation exceeds the total degrees of freedom present in the signals.

A special case of (1) where the constituents $\{\bm{X}_{k}\}_{k=1}^{s}$ are rank-one matrices has been analyzed from a theoretical and a numerical viewpoint in [14]. There, the authors investigate nuclear norm minimization (NNM) for the demixing problem under structured measurements of practical interest, and show that $m\gtrsim s^{2}n$ number of measurements are sufficient for NNM to reliably extract each constituent $\bm{X}_{k}$ from the single measurement vector $\bm{y}$. It is also worth noting that an improvement of the theoretical analysis in [14] has been announced in [21].

A limitation of the numerical methods proposed in [14] and [18] is that the resulting semidefinite program is computationally rather expensive to solve for medium-size and large-size problems. Some applications require numerical efficiency, in which case a different numerical approach is needed. The goal of this paper is to develop numerically efficient algorithms that can solve the nonlinear low rank matrix demixing problem without resorting to convex optimization, and meanwhile to provide competitive theoretical recovery guarantees for the proposed algorithms. Closest to the goal of this paper is arguably a very recent paper by Ling and one of the authors [15]. There, the authors consider the same setup as in [14], and propose a non-convex regularized gradient-descent based method. The differences to this paper are that (i) [15] is specialized to the joint blind deconvolution-demixing setting and is not designed for (1), where the unknown matrices $\bm{X}_{k}$ are general rank-$r$ matrices and the linear operators $\mathcal{A}_{k}$ are more general sensing matrices; (ii) while both algorithms fall in the realm of non-convex optimization, [15] uses a gradient descent method and this paper uses thresholding-based methods; (iii) the theoretical analysis in this paper does not apply to the joint blind deconvolution-demixing setting in [15].

1.2 Outline and Notation

The remainder of the paper is organized as follows. The numerical algorithms and their theoretical guarantees for the demixing problem are presented in Section 2. We also introduce an Amalgam form of restricted isometry property in Section 2 around which our theoretical analysis revolves. In Section 3 we test our algorithms on a variety of numerical examples. The proofs of the theoretical results are presented in Section 4. We conclude this paper with some potential future directions in Section 5.

Throughout the paper we use the following notational conventions. We denote vectors by bold lowercase letters and matrices by bold uppercase letters. In particular, we fix $\{\bm{X}_{k}\}_{k=1}^{s}$ as the target matrices and $\{\bm{X}_{k,l}\}_{k=1}^{s}$ as the iterates of the algorithms. For a matrix $\bm{Z}$, we use $\left\|\bm{Z}\right\|$ and $\left\|\bm{Z}\right\|_{F}$ to denote its spectral norm and Frobenius norm, respectively. For both vectors and matrices, $\bm{z}^{T}$ and $\bm{Z}^{T}$ denote their transpose while $\bm{z}^{*}$ and $\bm{Z}^{*}$ denote their conjugate transpose. The inner product of two matrices $\bm{Z}_{1}$ and $\bm{Z}_{2}$ is defined as $\left\langle\bm{Z}_{1},\bm{Z}_{2}\right\rangle=\operatorname{trace}(\bm{Z}_{1}^{*}\bm{Z}_{2})$.

2 Amalgam-RIP and Algorithms

2.1 Problem Setup and Amalgam-RIP

As outlined in the introduction, we want to solve the following demixing problem:

[TABLE]

where $\{\mathcal{A}_{k}\}_{k=1}^{s}$ is a set of linear operators mapping $n\times n$ matrices to $m$-dimensional vectors. Let $\{\bm{A}_{k,p}~{}|~{}1\leq k\leq s,1\leq p\leq m\}$ be a set of measurement matrices. We can write $\mathcal{A}_{k}(\bm{Z})$ explicitly as

[TABLE]

For conciseness, we only present our algorithms and theoretical results for $n\times n$ real matrices, but it is straightforward to modify them for $n_{1}\times n_{2}$ complex matrices.

We will propose two iterative hard thresholding algorithms for low rank matrix demixing. A question of central importance is how many measurements are needed so that the algorithms can successfully extract all the constituents $\{\bm{X}_{k}\}_{k=1}^{s}$ from $\bm{y}$. Since each $\bm{X}_{k}$ is determined by $r(2n-r)$ parameters [5], we need at least $m\geq sr(2n-r)$ many measurements for the problem to be well-posed. One of our goals is to show that the proposed algorithms can succeed already if the number of measurements is close to this information-theoretic minimum. As is common for hard thresholding algorithms in compressed sensing and low rank matrix recovery, the convergence analysis here will rely on some form of restricted isometry property (RIP). The notion most fitting to the demixing setting requires the RIP in an amalgamated form.

Definition 2.1** (Amalgam-RIP).**

The set of linear operators $\{\mathcal{A}_{k}\}_{k=1}^{s}$ satisfy the Amalgam-RIP (ARIP) with the parameter $\delta_{r}$ if

[TABLE]

holds for all the matrices $\bm{Z}_{k},~{}k=1,\cdots,s$ of rank at most $r$.

Fix an index $k$ and take $\bm{Z}_{k^{\prime}}=0$ for all $k^{\prime}\neq k$. The ARIP implies that

[TABLE]

hold for all matrices $\bm{Z}_{k}$ of rank at most $r$, which is indeed the RIP introduced in [19] for low rank matrix recovery. When the measurements matrices in (4) have i.i.d $\mathcal{N}(0,1)$ entries, we can show that $\{\mathcal{A}_{k}\}_{k=1}^{s}$ satisfy the ARIP with overwhelmingly high probability provided the number of measurements is proportional to the number of degrees of freedom within the constituents $\{\bm{X}_{k}\}_{k=1}^{s}$.

Theorem 2.1**.**

If $\{\bm{A}_{k,p}~{}|~{}1\leq k\leq s,1\leq p\leq m\}$ is a set of standard Gaussian matrices, then $\{\mathcal{A}_{k}\}_{k=1}^{s}$ satisfy the ARIP with the parameter $\delta$ with probability at least $1-2\exp(-m\delta^{2}/64)$ provided

[TABLE]

where $C>0$ is an absolute numerical constant.

2.2 Iterative Hard Thresholding Algorithms

In this subsection, we present two iterative hard thresholding algorithms for the low rank demixing problem. The algorithms are developed by targeting the following rank constraint problem directly

[TABLE]

The first algorithm, simply referred to as iterative hard thresholding (IHT), is presented in Algorithm 1. In each iteration, the algorithm first computes the residual and then updates each constituent separately by projected gradient descent. The search direction with respect to each constituent is computed as the negative gradient descent direction of the objective function when the other constituents are fixed. The hard thresholding operator $\mathcal{H}_{r}$ in the algorithm returns the best rank $r$ approximation of a matrix. The stepsize within the line search can be either fixed or computed adaptively. In Algorithm 1, we compute it as the steepest descent stepsize along the projected gradient descent direction $\mathcal{P}_{T_{k,l}}(\bm{G}_{k,l})$, where $T_{k,l}$ is the tangent space of the rank $r$ matrix manifold at the current estimate $\bm{X}_{k,l}$. Let $\bm{X}_{k,l}=\bm{U}_{k,l}\bm{\Sigma}_{k,l}\bm{V}_{k,l}^{*}$ be the singular value decomposition of $\bm{X}_{k,l}$. The tangent space $T_{k,l}$ consists of matrices which share the same column or row subspaces with $\bm{X}_{k,l}$ [1],

[TABLE]

Note that the inner loop within the algorithm is fully parallel since we update each constituent separately.

Theorem 2.2**.**

Assume the linear operators $\{\mathcal{A}_{k}\}_{k=1}^{s}$ satisfy the ARIP with the parameter $\delta_{3r}$. Define

[TABLE]

If $\gamma_{1}<1$, then the iterates of IHT satisfy

[TABLE]

*It follows immediately that the iterates of IHT converge linearly to the underlying constituents provided $\delta_{3r}<1/5$. *

If $\{\mathcal{A}_{k}\}_{k=1}^{s}$ consist of standard Gaussian measurement matrices, then Theorem 2.2 together with Theorem 2.1 implies that the number of necessary measurements for IHT to achieve successful recovery is $O(snr)$ which is optimal. In addition, we can establish the robustness of IHT against additive noise by considering the model $\bm{y}=\sum_{k=1}^{s}\mathcal{A}_{k}(\bm{X}_{k})+\bm{e}$, where $\bm{e}$ represents a noise term.

Theorem 2.3**.**

Assume $\bm{y}=\sum_{k=1}^{s}\mathcal{A}_{k}(\bm{X}_{k})+\bm{e}$ and the linear operators $\{\mathcal{A}_{k}\}_{k=1}^{s}$ satisfy the ARIP with parameter $\delta_{3r}$. If $\gamma_{1}<1$, then the iterates of IHT satisfy

[TABLE]

where

[TABLE]

The application of the hard thresholding operator in Algorithm 1 requires the singular value decomposition (SVD) of an $n\times n$ matrix in each iteration, which is computationally expensive for unstructured matrices. Inspired by the work of Riemannian optimization for low rank matrix reconstruction in [23, 27, 26], we propose to accelerate IHT by updating each constituent along a projected gradient descent direction, see fast iterative hard thresholding (FIHT) described in Algorithm 2. The key difference between Algorithms 1 and 2 lies in Step 3. Instead of updating $\bm{X}_{k,l}$ along the gradient descent direction $\bm{G}_{k,l}$ as in IHT, FIHT updates $\bm{X}_{k,l}$ along the projected gradient descent direction $\mathcal{P}_{T_{k,l}}(\bm{G}_{k,l})$, where $T_{k,l}$ is the tangent space defined in (6), and then followed by the hard thresholding operation.

Let $\bm{W}_{k,l}=\bm{X}_{k,l}+\alpha_{k,l}\mathcal{P}_{T_{k,l}}(\bm{G}_{k,l})$. One can easily observe that $\bm{W}_{k,l}\in T_{k,l}$ and $\operatorname{rank}(\bm{W}_{k,l})\leq 2r$. Moreover, $\bm{W}_{k,l}$ obeys the following decomposition [27]:

[TABLE]

where $\bm{Q}_{1}$ and $\bm{Q}_{2}$ are two $n\times r$ orthonormal matrices such that $\bm{Q}_{1}\perp\bm{U}_{k,l}$ and $\bm{Q}_{2}\perp\bm{V}_{k,l}$, and $\bm{M}_{k,l}$ is a $2r\times 2r$ matrix of the form

[TABLE]

Consequently, the SVD of $\bm{W}_{k,l}$ can be obtained from the SVD of the $2r\times 2r$ matrix $\bm{M}_{k,l}$. Therefore, the introduction of the extra projection $\mathcal{P}_{T_{k,l}}$ for the search direction can reduce the computational complexity of the partial SVD of $\bm{W}_{k,l}$ from $O(n^{2}r)$ flops to $O(nr^{2}+r^{3})$ flops. The recovery guarantee of FIHT can also be established in terms of the ARIP.

Theorem 2.4**.**

*Assume the linear operators $\{\mathcal{A}_{k}\}_{k=1}^{s}$ satisfy the ARIP with the parameter $\delta_{3r}$. Define *

[TABLE]

where $\operatorname{\sigma_{\max}}:=\max_{k}\operatorname{\sigma_{\max}}(\bm{X}_{k})$ and $\operatorname{\sigma_{\min}}:=\max_{k}\operatorname{\sigma_{\min}}(\bm{X}_{k})$. If $\gamma_{2}<1$, then the iterates of Algorithm 2 satisfy

[TABLE]

It follows that the iterates of FIHT converge linearly to the underlying constituents provided

[TABLE]

Assume $\{\mathcal{A}_{k}\}_{k=1}^{s}$ consist of standard Gaussian matrices. In contrast to the $O(srn)$ necessary measurements for IHT to achieve successful recovery, Theorem 2.4 implies that FIHT requires $O((\operatorname{\sigma_{\max}}/\operatorname{\sigma_{\min}})^{2}s^{2}nr^{2})$ measurements which is not optimal. However, numerical simulations in Section 3 suggest that FIHT needs fewer measurements than IHT to be able to successfully reconstruct $s$ rank-$r$ matrices.

3 Numerical Simulations and Applications

We evaluate the empirical performance of IHT and FIHT on simulated problems as well as application examples from quantum states demxing and the Internet of Things. In the implementation IHT and FIHT are terminated when a maximum of $500$ iterations is met or the relative residual is small,

[TABLE]

They are considered to have successfully recovered a set of matrices $\{\bm{X}_{k}\}_{k=1}^{s}$ if the returned approximations $\{\bm{X}_{k,l}\}_{k=1}^{s}$ satisfy

[TABLE]

All the random tests are repeated ten times in this section.

3.1 Phase Transitions under Gaussian Measurements

The ARIP-based theoretical recovery guarantees in Theorems 2.2 and 2.4 are worst-case analysis which are uniform over all rank $r$ matrices. However, these conditions are highly pessimistic when compared with average-case empirical observations. In this subsection we investigate the empirical recovery performance of IHT and FIHT under Gaussian measurements. The tests are conducted on $n\times n$ rank-$r$ matrices with $n=50$ and $r=5$, and the number of constituents $s$ varies from $1$ to $7$. For a fixed value of constituents, a number of equispaced values of $m$ are tested.

We first test IHT and FIHT on well-conditioned matrices which are formed via $\bm{X}_{k}=\bm{L}_{k}\bm{R}_{k}^{*}$, where $\bm{L}_{k}\in\mathbb{R}^{n\times r}$, $\bm{R}_{k}\in\mathbb{R}^{n\times r}$, and their entries are drawn i.i.d. from the standard normal distribution. The phase transitions of IHT and FIHT on well-conditioned matrices are presented in Figures 1a and 1b, respectively. Figure 1a shows a linear correlation between the number of measurements $m$ and the number of constituents $s$ for the successful recovery of IHT. Though the recovery guarantee of FIHT in Theorem 2.4 is worse than that of IHT in Theorem 2.2, Figure 1b shows that FIHT requires fewer measurements than IHT to achieve successful recovery with high probability. In particular, $m\approx 3.3sr(2n-r)$ is sufficient for FIHT to successfully reconstruct a set of $s$ low rank matrices for all the tested values of $s$.

To explore how the condition number of the test matrices impact the recovery of FIHT, we further test FIHT on ill-conditioned matrices. The matrices are computed as $\bm{X}_{k}=\bm{U}_{k}\bm{\Sigma}\bm{V}_{k}^{*}$, where $\bm{U}_{k}$ and $\bm{V}_{k}$ are two $50$ by $5$ random orthonormal matrices, and $\bm{\Sigma}$ is a $5$ by $5$ diagonal matrix with the diagonal entries being $5$ uniformly spaced values from $1$ to $1000$. So the condition number of the test matrices is $1000$. The phase transition plot of FIHT on ill-conditioned matrices is presented in Figure 1c. The figure shows that FIHT requires more measurements to reconstruct the test ill-conditioned matrices than to reconstruct the well-conditioned matrices when $s\geq 5$. However, FIHT is still able to reconstruct all the test ill-conditioned matrices when $m\approx 3.8sr(2n-r)$. Hence, there is no strong evidence from Figure 1c to support that the number of measurements for the successful recovery of FIHT relies on the square of the condition number of the target matrices, which may suggest the possibility of improving the result in Theorem 2.4.

In the reminder of this section, we focus on the empirical performance of FIHT for its superiority and flexibility. Though we have only established the recovery guarantee of FIHT for real matrices under the Gaussian measurements, the algorithm is equally effective for complex matrices and other types of measurements. In the next two subsections, we study the applications of FIHT in quantum states demixing and the Internet of Things.

3.2 An Application in Quantum States Demixing

As stated in the introduction, the state of a quantum system can be represented by a low rank, positive semidefinite density matrix of unit trace. Next, we present a stylized application of FIHT to the demixing of quantum states. When the target constituents are $n$ by $n$ rank-$r$ positive semidefinite matrices of unit trace, FIHT in Algorithm 2 can be modified accordingly to preserve the matrix structures, see Algorithm 3 for complete description.

Let $\bm{X}_{k,l}$ be the current rank-$r$ estimate which is Hermitian, positive semidefinite, and so admits an eigenvalue decomposition (equivalent to its SVD) $\bm{X}_{k,l}=\bm{U}_{k,l}\bm{\Lambda}_{k,l}\bm{U}_{k,l}^{*}$ with non-negative eigenvalues. There are two key differences between Algorithms 2 and 3. Firstly, the tangent space can be defined as

[TABLE]

so that all the matrices in $T_{k,l}$ are Hermitian. Therefore after $\bm{X}_{k,l}$ is updated along the projected gradient descent direction, $\bm{W}_{k,l}=\bm{X}_{k,l}+\alpha_{k,l}\mathcal{P}_{T_{k,l}}(\bm{G}_{k,l})$ remains Hermitian. In addition, the eigenvalue decomposition of $\bm{W}_{k,l}$ can also be computed using $O(nr^{2})$ flops as in the computation of the SVD of the non-symmetric matrices in Algorithm 2, since $\bm{W}_{k,l}$ has the following decomposition:

[TABLE]

where $\begin{bmatrix}\bm{U}_{k,l}&\bm{Q}\end{bmatrix}$ is an $n\times 2r$ orthonormal matrix and $\bm{M}_{k,l}$ is a $2r\times 2r$ Hermitian matrix of the form

[TABLE]

Secondly, the hard thresholding operator $\mathcal{H}_{r}$ is replaced by the projection $\mathcal{P}_{\Pi}$, followed by another projection $\mathcal{P}_{\Delta}$, where $\mathcal{P}_{\Pi}$ projects $\bm{W}_{k,l}$ onto the set of low rank positive semidefinite matrices

[TABLE]

and $\mathcal{P}_{S}$ projects the eigenvalues of $\mathcal{P}_{\Pi}(\bm{W}_{k,l})$ onto the simplex

[TABLE]

Notice that the inertia of a Hermitian matrix $\bm{H}$, denoted by $\mbox{inertia}(\bm{H})$, is a triple of integer numbers ($i_{p},i_{n},i_{z}$), where $i_{p}$, $i_{n}$, $i_{z}$ are the number of positive, negative, and zero eigenvalues of $\bm{H}$. From [9, Lemma A.15], we know that $\mbox{inertia}(\bm{M}_{k,l})=\mbox{inertia}(\bm{N}^{*}\bm{M}_{k,l}\bm{N})+(t,t,r-t)$, where $t=\operatorname{rank}(\bm{R})$ and $\bm{N}$ is an $r\times(r-t)$ matrix whose columns forms a null space of $\bm{R}$. This implies that $\bm{M}_{k,l}$ has at most $r$ positive eigenvalues, so does $\bm{W}_{k,l}$. Therefore by [11, Theorem 2.1], we can compute $\mathcal{P}_{\Pi}(\bm{W}_{k,l})$ by only keeping the positive eigenvalues and the corresponding eigenvectors of $\bm{M}_{k,l}$. To compute $\mathcal{P}_{S}(\cdot)$, we use an algorithm proposed in [7], in which the computational cost is dominated by sorting a vector of length $r$.

We evaluate the recovery performance of FIHT for quantum states demixing under the Pauli measurements, where the measurement matrices in (4) are tensor products of the Pauli matrices. That is,

[TABLE]

where

[TABLE]

The target matrices are computed as $\bm{X}_{k,l}=\frac{1}{r}\bm{U}_{k,l}\bm{U}_{k,l}^{*}$, where $\bm{U}_{k,l}$ are random $n\times r$ complex orthonormal matrices. We conduct the experiments with $q=6$, $r=1\mbox{ and }3$, and $a$ is selected uniformly at random. When $r=1$, it corresponds to demixing the pure states in quantum mechanics. The phase transitions of FIHT under the Pauli measurements is presented in Figure 2, which shows that FIHT is equally effective for quantum states demixing and $m\approx 3.5sr(2n-r)$ is sufficient for successful demixing all the tested values of $s$ from $1$ to $7$.

3.3 An Application from the Internet-of-Things

Demixing problems of the form (1) are also expected to play an important role in the future IoT. In this subsection we conduct some numerical simulations for a demixing problem related to the IoT. Before we proceed to the simulations we describe in more detail how such a demixing problem does arise in the IoT and how it is connected to the demixing of low rank matrices.

In mathematical terms we are dealing—in somewhat simplified form—with the following problem. We consider a scenario where many different users/devices communicate with a base station, see Figure 3 for an illustration of the described setup. Let $\bm{x}_{k}\in\mathbb{C}^{n_{1}}$ be the data that are supposed to be transmitted by the $k$-th user. Before transmission we encode $\bm{x}_{k}$ by computing $\bm{z}_{k}=\bm{C}_{k}\bm{x}_{k}$, where $\bm{C}_{k}\in\mathbb{C}^{m\times n_{1}}$ is a fixed precoding matrix of the $k$-th user. The encoding matrice $\bm{C}_{k}$ differs from user to user but is known to the receiver. Let $\bm{g}_{k}$ be the associated impulse response of the communication channel between the $k$-th user and the base station. Since we assume that the channel does not change during the transmission of the signal $\bm{z}_{k}$, $\bm{g}_{k}$ acts as convolution operator, i.e., the signal arriving at the base station becomes $\bm{z}_{k}\ast\bm{g}_{k}$ (ignoring additive noise for simplicity). However, the channels $\bm{g}_{k}$ are not known to the base station (or to the transmitter), since this is the overhead information that we want to avoid transmitting, as it may change from transmission to transmission. The signal recorded at the base station is given by a mixture of all transmitted signals convolved with their respective channel impulse responses, i.e.,

[TABLE]

The goal is to demix the received signal $\bm{y}$ and recover the signals $\bm{x}_{k}$ of the individual users.

But how is (7) related to our setup in (1)? We first note that in all wireless communication scenarios of relevance the impulse response $\bm{g}_{k}$ is compactly supported. In engineering terms the size of the support of $\bm{g}_{k}$, denoted here by $n_{2}$, is called its maximum delay spread. In practice, the length of the transmitted signals $m$ is typically (much) larger than the maximum delay spread $n_{2}$. Thus, $\bm{g}_{k}$ is a $m$-dimensional vector which is formed by the channel impulse responses being padded with zeroes. We can write $\bm{g}_{k}$ as $\bm{g}_{k}=[\bm{h}_{k}^{T},0,\dots,0]^{T}$, where $\bm{h}_{k}$ is the non-zero padded impulse response of length $n_{2}$. Let $\bm{F}$ denote the $m\times m$ unitary Discrete Fourier Transform matrix and let $\bm{B}$ be the $m\times n_{2}$ matrix consisting of the first $n_{2}$ columns of $\bm{F}$. By the circular convolution theorem, we can express the convolution between $\bm{z}_{k}$ and $\bm{g}_{k}$ in the Fourier domain via

[TABLE]

where $\odot$ denotes the componentwise product. Let $\bm{c}_{k,p}^{*}$ represent the $p$-th row of $\bm{F}\bm{C}_{k}$ and let $\bm{b}_{p}^{*}$ represent the $p$-th row of $\bm{B}$. A simple calculation shows that [14]

[TABLE]

where $\bm{X}_{k}:=\bm{x}_{k}\bar{\bm{h}}_{k}^{*}$ is an $n_{1}\times n_{2}$ rank-one matrix containing the unknown signal and the unknown channel of the $k$-th user, while $\bm{c}_{k,p}$ and $\bm{b}_{p}^{\ast}$ are known vectors (since $\bm{C}_{k}$ and $\bm{B}$ are known). Hence, at the base station the observed signal $\bm{y}=\sum_{k=1}^{s}\bm{z}_{k}\ast\bm{g}_{k}$ can be expressed in the Fourier domain as $\frac{1}{\sqrt{m}}\hat{\bm{y}}=\sum_{k=1}^{s}\mathcal{A}_{k}(\bm{X}_{k})$, which is of the same form as (1) (modulo replacing $\frac{1}{\sqrt{m}}\hat{\bm{y}}$ with $\bm{y}$). The measurement matrices in this scenario are given by $\bm{A}_{k,p}=\bm{c}_{k,p}\bar{\bm{b}}_{p}^{*}$.

We test FIHT (Algorithm 2) for the demixing problem described in (8) with $n_{1}=128$ and $n_{2}=32$, which are of interest in practice. As a result, the constituents are non-square matrices. Two types of encoding matrices are tested: a) the entries of $\bm{C}_{k}$ are standard complex Gaussian variables; and b) $\bm{C}_{k}=\bm{D}_{k}\bm{H}$, where $\bm{H}$ is the first $n_{2}$ columns of an $m\times m$ Hadamard matrix, and then premultiplied by an $m\times m$ diagonal random sign matrix $\bm{D}$. The target signals $\bm{x}_{k}$ and $\bm{h}_{k}$ are complex Gaussian random vectors. The phase transition plot for $\bm{C}_{k}$ being a Gaussian matrix is presented in Figure 4a, where the number of constituents varies from $1$ to $7$. Figure 4a displays a linear correlation between the number of necessary measurements and the number of the constituents, and it shows that $m\approx 3s(n_{1}+n_{2})$ number of measurements are sufficient for successful demixing of the transmission signals and the impulse responses. When $\bm{C}_{k}$ is formed from a partial Hadamard matrix, we only test $m=2^{10},2^{11},~{}2^{12}$, and $2^{13}$ due to the Hadamard conjecture [10], but a larger range of $s$ are tested. The phase transition plot presented in Figure 4b shows that FIHT is still able to successfully demix the transmission signals and the impulse responses when $m\approx 3s(n_{1}+n_{2})$.

3.4 Robustness to Additive Noise

We explore the robustness of FIHT against additive noise by conducting tests on the problem instances set up in the last three subsections. We have established in Theorem 2.3 that Algorithm 1 (IHT) is stable in the presence of additive noise, but could not provide a similar theoretical robustness guarantee for FIHT because of the more involved nature of the current convergence analysis framework of FIHT. Yet, empirical evidence clearly indicates that FIHT also shows the desirable stability vis-à-vis noise, as we will demonstrate in this subsection.

We assume that the measurement vectors $\bm{y}$ in the tests are corrupted by

[TABLE]

where $\bm{w}$ is an $m\times 1$ standard Gaussian vector, either real or complex up to the testing environment, and $\sigma$ takes nine different values of $\sigma$ from $10^{-4}$ to $1$. For each problem instance, two values of the number of measurements are tested. The average relative reconstruction error out in dB against the signal to noise ratio (SNR) is plotted in Figure 5. The figure shows that the relative reconstruction error scales linearly with the noise levels under different measurement schemes. Moreover, as desired, the relative reconstruction error decreases linearly on a log-log scale as the number of measurements increases.

3.5 A Rank Increasing Heuristic

While in some applications the rank of the target matrices is known a priori, for example in pure states demxing and the IoT application, we may not know the matrix rank exactly in other applications. Here, we suggest a rank increasing heuristic when the rank of the matrices is not known. Starting from rank one, we execute FIHT until the residual is not decreasing significantly; and then we successively increase the rank and repeat this procedure until the residual is sufficiently small. We test this heuristic using problems instances from Subsections 3.1 and 3.2. The relative residual plotted against the number of iterations is presented in Figure 6. The jumps in the error curves correspond to the iteration steps when the algorithm increases the assumed rank of the unknown matrices by 1.

4 Proofs

4.1 Proof of Theorem 2.1

Proof.

Let $\{\bm{Z}_{k}\}_{k=1}^{s}$ be a set of fixed matrices with $\operatorname{rank}(\bm{Z}_{k})\leq r$ for all $1\leq k\leq s$. One can easily see that

[TABLE]

Since for any $1\leq p\leq m$, $\bm{A}_{1,p},\bm{A}_{2,p},\cdots,\bm{A}_{s,p}$ are independent from each other, one has

[TABLE]

Therefore,

[TABLE]

where $\xi_{p}^{2},~{}p=1,\cdots,m$ are i.i.d $\chi^{2}$ variables. Consequently,

[TABLE]

and the application of the Bernstein inequality for Chi-squared variables [24] gives

[TABLE]

Taking $t=\delta\cdot m\sum_{k=1}^{s}\left\|\bm{Z}_{k}\right\|_{F}^{2}$ gives

[TABLE]

for $0<\delta<1$. This implies for fixed $\left\{\bm{Z}_{k}\right\}_{k=1}^{s}$, we have that

[TABLE]

holds with probability at least $1-2e^{-m\delta^{2}/8}$.

The rest of the proof follows essentially from the proof of Theorem 2.3 in [5]. Due to the homogeneity of the problem with respect to $\sum_{k=1}^{s}\left\|\bm{Z}_{k}\right\|_{F}^{2}$, we only need to show (10) for the set of matrices in

[TABLE]

Based on the SVD of each $\bm{Z}_{k}$, we have

[TABLE]

Define $\mathcal{Q}_{n,r}=\{\bm{Q}\in\mathbb{R}^{n\times r}:~{}\bm{Q}^{*}\bm{Q}=\bm{I}\}$. It has been shown in [5] that $\mathcal{Q}_{n,r}$ has a $\epsilon$-net $\overline{\mathcal{Q}}_{n,r}\subset\mathcal{Q}_{n,r}$ of cardinality $|\overline{\mathcal{Q}}_{n,r}|\leq(3/\epsilon)^{nr}$ under the $\left\|\cdot\right\|_{1,2}$ norm. In other words, for any $\bm{Q}\in\mathcal{Q}_{n,r}$, there exists a $\overline{\bm{Q}}\in\overline{\mathcal{Q}}_{n,r}$ such that

[TABLE]

Let $\mathcal{D}$ be the set of $rs\times rs$ diagonal matrices with unit norm diagonal entries. It is well known that $\mathcal{D}$ has a $\epsilon$-net $\overline{\mathcal{D}}\subset\mathcal{D}$ of cardinality $|\overline{\mathcal{D}}|\leq(3/\epsilon)^{rs}$ [24]. Define

[TABLE]

We are going to show that $\overline{\mathcal{S}}_{r}$ is a $3\epsilon$-net of $\mathcal{S}_{r}$ with cardinality $|\overline{\mathcal{S}}_{r}|\leq(3/\epsilon)^{(2n+1)rs}$.

First note that for any $[\bm{U}_{1}\bm{\Sigma}_{1}\bm{V}_{1}^{*},\cdots,\bm{U}_{s}\bm{\Sigma}_{s}\bm{V}_{s}^{*}]^{T}$, there exists a $[\overline{\bm{U}}_{1}\overline{\bm{\Sigma}}_{1}\overline{\bm{V}}_{1}^{*},\cdots,\overline{\bm{U}}_{s}\overline{\bm{\Sigma}}_{s}\overline{\bm{V}}_{s}^{*}]^{T}\in\overline{\mathcal{S}}_{r}$ such that $\left\|\bm{U}_{k}-\overline{\bm{U}}_{k}\right\|_{1,2}\leq\epsilon$, $\left\|\bm{V}_{k}-\overline{\bm{V}}_{k}\right\|_{1,2}\leq\epsilon$, and $\sqrt{\sum_{k=1}^{s}\left\|\bm{\Sigma}_{k}-\overline{\bm{\Sigma}}_{k}\right\|_{F}^{2}}\leq\epsilon$. Thus,

[TABLE]

Since

[TABLE]

we have $\operatorname{{\textbf{\mbox{I}}}}_{1}\leq\epsilon$. Similarly, the third term can be bounded as $\operatorname{{\textbf{\mbox{I}}}}_{3}\leq\epsilon$. To bound $\operatorname{{\textbf{\mbox{I}}}}_{2}$, note that

[TABLE]

So we also have $\operatorname{{\textbf{\mbox{I}}}}_{2}\leq\epsilon$. Combining the bounds for $\operatorname{{\textbf{\mbox{I}}}}_{1},~{}\operatorname{{\textbf{\mbox{I}}}}_{2}$ and $\operatorname{{\textbf{\mbox{I}}}}_{3}$ together implies that $\overline{\mathcal{S}}_{r}$ is a $3\epsilon$-net of $\mathcal{S}_{r}$ with cardinality $|\overline{\mathcal{S}}_{r}|\leq(3/\epsilon)^{(2n+1)rs}$. Therefore,

[TABLE]

provided

[TABLE]

So for all $\left\{\overline{\bm{Z}}_{k}\right\}\in\overline{\mathcal{S}}_{k}$, we have

[TABLE]

with probability at least $1-2e^{-m\delta^{2}/64}$.

Define

[TABLE]

In the following, we will take $\epsilon=\delta/(12\sqrt{2})$. So $\overline{\mathcal{S}}_{r}$ is a $\delta/(4\sqrt{2})$-net of $\mathcal{S}_{r}$; and for any $\{\bm{Z}_{k}\}_{k=1}^{s}\in\mathcal{S}_{r}$, there exist $\{\overline{\bm{Z}}_{k}\}_{k=1}^{s}\in\overline{\mathcal{S}}_{r}$ such that

[TABLE]

Thus,

[TABLE]

Since $\bm{Z}_{k}-\overline{\bm{Z}}_{k}$ is a matrix of rank at most $2r$, we can decompose it as

[TABLE]

where $\left\langle\bm{Y}_{k}^{1},\bm{Y}_{k}^{2}\right\rangle=0$, $\operatorname{rank}(\bm{Y}_{k}^{i})\leq r~{}(i=1,2)$, and $\left\|\bm{Y}_{k}^{1}\right\|_{F}^{2}+\left\|\bm{Y}_{k}^{2}\right\|_{F}^{2}=\left\|\bm{Z}_{k}-\overline{\bm{Z}}_{k}\right\|_{F}^{2}$. Then it follows that

[TABLE]

where the second line follows from $\left\{\bm{Y}_{k}^{i}/\sqrt{\sum_{k=1}^{s}\left\|\bm{Y}_{k}^{i}\right\|_{F}^{2}}\right\}_{k=1}^{s}\in\mathcal{S}_{r}$ for $i=1,2$. Inserting (12) into (11) gives

[TABLE]

which implies $\kappa_{r}\leq 1+\delta$ and

[TABLE]

for all $\{\bm{Z}_{k}\}_{k=1}^{s}\in{\mathcal{S}_{r}}$. The lower bound can be obtained in the following way:

[TABLE]

where the second line follows from (12) and the bound for $\kappa_{r}$. The upper and lower bounds in (13) and (14) can be easily transferred into the desired ARIP bounds by a change of variables. ∎

4.2 Proof of Theorem 2.2

Proof.

Let $\bm{U}_{k,l+1}$, $\bm{U}_{k,l}$, and $\bm{U}_{k}$ be the left singular vectors of $\bm{X}_{k,l+1}$, $\bm{X}_{k,l}$ and $\bm{X}_{k}$, respectively. Let $\bm{Q}_{k,l}$ be an $n\times 3r$ orthonormal matrix which spans the union of the column spaces of $\bm{U}_{k,l+1}$, $\bm{U}_{k,l}$, and $\bm{U}_{k}$. We have $\bm{X}_{k,l+1}=\mathcal{P}_{\bm{U}_{k,l+1}}(\bm{X}_{k,l}+\alpha_{k,l}\bm{G}_{k,l})$. Moreover, by comparing the following two equalities,

[TABLE]

we have

[TABLE]

since $\bm{X}_{k,l+1}$ is the best rank $r$ approximation of $\bm{X}_{k,l}+\alpha_{k,l}\bm{G}_{k,l}$. Thus,

[TABLE]

Consequently we have

[TABLE]

where in the third line, we have used the fact $\mathcal{P}_{\bm{Q}_{k,l}}(\bm{X}_{k,l}-\bm{X}_{k})=(\bm{X}_{k,l}-\bm{X}_{k})$ for all $1\leq k\leq s$. In order to bound $\operatorname{{\textbf{\mbox{I}}}}_{1}$, we first rewrite it into the following matrix-vector product form

[TABLE]

where

[TABLE]

So it suffices to bound the spectral norm of $\mathcal{A}_{\alpha}$. Denote $\mathcal{A}_{\alpha}$ by $\mathcal{A}$ when $\alpha_{k,l}=1$ for all $1\leq k\leq s$ in $\mathcal{A}_{\alpha}$. Let $\bm{Y}=[\bm{Y}_{1},\cdots,\bm{Y}_{s}]^{T}$. Since $\mathcal{A}$ is self-adjoint, we have

[TABLE]

where the third line follows from the ARIP and the fact $\mathcal{P}_{\bm{Q}_{k,l}}(\bm{Y}_{k}),~{}k=1,\cdots,s$ are matrices of rank at most $3r$. It follows that

[TABLE]

where in the third line we use the following ARIP bound for $\alpha_{k,l}$

[TABLE]

Combining the spectral of $\mathcal{A}_{\alpha}$ together with (15) and (16) gives

[TABLE]

which completes the proof of the theorem. ∎

4.3 Proof of Theorem 2.3

Proof.

To prove the convergence of IHT under additive noise, we only need to modify (15) slightly as follows:

[TABLE]

The proof of Theorem 2.2 shows that

[TABLE]

while $\operatorname{{\textbf{\mbox{I}}}}^{\prime}_{12}$ can be bounded as

[TABLE]

where in the third line we utilize the bound for $\alpha_{k,l}$ in (17) and in the last line we utilize the definition of the ARIP. Therefore,

[TABLE]

The proof is complete after we apply the above inequality recursively. ∎

4.4 Proof of Theorem 2.4

Proof.

Let $\bm{W}_{k,l}=\bm{X}_{k,l}+\alpha_{k,l}\mathcal{P}_{T_{k,l}}(\bm{G}_{k,l})$. The following inequality holds:

[TABLE]

where the inequality follows from

[TABLE]

This leads to

[TABLE]

Following the same argument for the bound of $\operatorname{{\textbf{\mbox{I}}}}_{1}$ in (15), we can bound $\operatorname{{\textbf{\mbox{I}}}}_{2}$ as

[TABLE]

by noting that all the matrices in $T_{k,l}$ are of rank at most $2r$.

Next, $\operatorname{{\textbf{\mbox{I}}}}_{3}$ can be bounded in the following way:

[TABLE]

where the third line follows from Lemma 5.1 in the appendix, and the fourth line follows from the ARIP bound for $\alpha_{k,l}$ in (17).

To bound $\operatorname{{\textbf{\mbox{I}}}}_{4}$, first note that the application of Lemma 5.2 in the appendix gives

[TABLE]

where $\operatorname{\sigma_{\min}}:=\min_{k}\operatorname{\sigma_{\min}}(\bm{X}_{k})$. Thus,

[TABLE]

Substituting the bounds for $\operatorname{{\textbf{\mbox{I}}}}_{2}$, $\operatorname{{\textbf{\mbox{I}}}}_{3}$ and $\operatorname{{\textbf{\mbox{I}}}}_{4}$ into (18) gives

[TABLE]

Assume

[TABLE]

By induction, we have

[TABLE]

Note that the initial guess is obtained by one-step hard thresholding with $\alpha_{k,l}=1$. So following the proof of Theorem 2.2 we have

[TABLE]

Therefore, one has $\rho<1$ if

[TABLE]

Since

[TABLE]

(19) holds if

[TABLE]

where $\operatorname{\sigma_{\max}}:=\max_{k}\operatorname{\sigma_{\max}}(\bm{X}_{k})$. ∎

5 Conclusion and Future Directions

We have presented the first computationally efficient algorithms that can extract low rank matrices from a sum of linear measurements. These algorithms have potential applications in areas such as wireless communications and quantum tomography. Numerical simulations show an empirical performance that is quite close to the information theoretic limit in terms of demixing ability with respect to the number of measurements.

At the same time, there are still a number of open questions for further directions. Firstly, the robustness analysis of FIHT needs to be addressed in the future. Secondly, our theoretical framework so far is still a bit restrictive, since it only yields close-to-optimal results for Gaussian measurement matrices, which are however of limited use in applications. Thus, one future challenge consist in establishing good Amalgam-RIP bounds for structured measurement matrices. Thirdly, it is also interesting to see whether the Amalgam-RIP can be used to analyse other approaches for low rank matrix demixing. In particular, we want to investigate whether the Amalgam form of RIP in (5), but restricted onto a local subspace, is sufficient for the guarantee analysis of the nuclear norm minimization studied in [14]. Finally, in this paper the Amalgam-RIP is defined for homogeneous data, i.e., matrices of low rank. It is likely that similar Amalgam-RIP can be established for heterogeneous data and then be used in the analysis of different reconstruction programs.

Appendix

Lemma 5.1**.**

Suppose $\left\langle\bm{Y}_{k},\bm{Z}_{k}\right\rangle=0$ and $\operatorname{rank}(\bm{Y}_{k}+\bm{Z}_{k})\leq c$ for all $1\leq k\leq s$. Then

[TABLE]

Proof.

Due to the homogeneity, we can assume $\sqrt{\sum_{k=1}^{s}\left\|\bm{Y}_{k}\right\|_{F}^{2}}=1,~{}\sqrt{\sum_{k=1}^{s}\left\|\bm{Z}_{k}\right\|_{F}^{2}}=1$. Since $\left\langle\bm{Y}_{k},\bm{Z}_{k}\right\rangle=0$, we have $\sum_{k=1}^{s}\left\|\bm{Y}_{k}+\bm{Z}_{k}\right\|_{F}^{2}=2$. So

[TABLE]

following from Def. 2.1. Then the application of the parallelogram identity implies

[TABLE]

which concludes the proof. ∎

Lemma 5.2** ([27, Lem. 4.1]).**

Let $\bm{X}_{k,l}=\bm{U}_{k,l}\bm{\Sigma}_{k,l}\bm{V}^{*}_{k,l}$ be rank $r$ matrix and $T_{k,l}$ be the tangent space of the rank $r$ matrix manifold at $\bm{X}_{k,l}$. Let $\bm{X}_{k}$ be another rank $r$ matrix. One has

[TABLE]