Measurement Manipulation of the Matrix Sensing Problem to Improve Optimization Landscape

iii) I also think that the presence of the term $\sqrt{n^2 / m}$ in their bound is very limiting. Does this not necessitate the number of measurement to be proportional to $n^2$, i.e. grow as the ambient dimension? Much of the compressed sensing literature is about avoiding this situation.

Thank you for bringing up an important issue. We could directly get the following result from the proof of Theorem 4:

$$\alpha||\operatorname{vec}{(M)}||\_2^2=\frac{1-\delta\_{s}}{\sigma\_1^2(\mathbf{A})}||\operatorname{vec}{(M)}||\_2^2\leq||P\mathbf{A}\operatorname{vec}{(M)}||\_2^2\leq \frac{1+\delta\_{s}}{\sigma\_m^2(\mathbf{A})}||\operatorname{vec}{(M)}||\_2^2 =\beta||\operatorname{vec}{(M)}||\_2^2$$

Since rescaling will not affect the landscape of the matrix sensing problem, we multiply all of the above inequalities by $\frac{2}{\alpha+\beta}$ and the tightest RIP bound we could have is $\frac{\beta-\alpha}{\beta+\alpha}$, which is

$$\frac{(1+\delta\_s)\sigma_1^2(\mathbf{A}) - (1-\delta\_s)\sigma\_m^2(\mathbf{A})}{(1+\delta\_s)\sigma\_1^2(\mathbf{A}) + (1-\delta\_s)\sigma_m^2(\mathbf{A})} = 1-\frac{2(1-\delta\_s)\sigma\_m^2(\mathbf{A})}{(1+\delta\_s)\sigma\_1^2(\mathbf{A}) + (1-\delta\_s)\sigma_m^2(\mathbf{A})}$$

If we plug in the bound in Assumption 1 that

\operatorname{Pr}\left\\{\sqrt{n^2/ m}-1-t \leq \sigma\_i(\mathbf{A}) \leq 1+\sqrt{n^2 / m}+t, i \in[m]\right\\} \geq 1-2 \exp \left(-mt^2 / 2\right), \quad \forall \epsilon>0.

For $0<t<\sqrt{\frac{n^2}{m}}-1$, $t$ can go to infinity while $\frac{t}{\sqrt{n^2/m}}\to 0$, we could have

$$\delta\_{s}(\mathcal{\tilde{A}}) \leq 1-\frac{2}{1+\frac{(1+\delta\_s)\sigma_1^2(\mathbf{A})}{ (1-\delta\_s)\sigma_m^2(\mathbf{A})}} \leq 1-\frac{2}{1+\frac{(1+\delta\_s)[\sqrt{n^2/ m}+t+1 ]^2}{ (1-\delta\_s)[\sqrt{n^2/ m}-t-1 ]^2}}$$

As $m,n\to \infty$, in the order of $m \gtrsim ns\log n$, we could have $\frac{[\sqrt{n^2/ m}+t+1 ]^2}{ [\sqrt{n^2/ m}-t-1 ]^2}\to 1$, and the right-hand side upper bound goes to $\delta\_s$. Thanks a lot for mentioning this concern, we will update the new bound in our modified paper. There is no need for $m$ to be on the order of $n^2$ and it should be $n\log n$.

3. Small inaccuricies

While the paper generally is well written, it would benefit from another round of proofreading. There are some inconsistencies in the notation: The norm $||A||\_\infty$ is never defined, on page 9, there is a term A(M) that should be a A(X), and so forth.

Thanks a lot for pointing that out. We will fix the typos and proofread the paper one more time.

2. Unclear interpretation of Theorems 4 and 5

i) From my point of view, Theorems 4 and 5 do little to explain the empirical success of the preconditioning strategy. As the authors point out, Theorem 4 only provides a better bound on the RIP constant if $\sigma_1(A)^2 \leq 1+\delta_s$. Since $\delta_s \geq \max \left(1-\sigma_m(A)^2, \sigma_1(A)^2-1\right)$, this can essentially only be the case if the singular values of $A$ are biased downward.

For the simulation experiments, we could see that, for uniform, correlated normal and poisson distribution, the original sensing operator has a RIP constant close to 1, where the Gaussian assumption violated. These are the main situations we care about, and we could see that the pre-conditioning has a great improvement on those distributions. The method will have improvement on those worst cases, and for the violation of $\sigma_1(A)^2 \leq 1+\delta_s$, the original $\delta_s$ is not too large to substantially harm the landscape.

ii) In particular, in the Gaussian setting that they consider in Theorem 5 , it is not. Indeed, $1-(1-\delta) /\left(1+\sqrt{n^2 / m}(1+\epsilon) \geq \delta\right.$ for all values of $\delta$-- i.e., the 'direct bound' that is given by previous results is better. In their proof, they seem to argue that since $A$ has the $\delta$-RIP, it also has the $2 \delta$-RIP, and then compare their bound to $2 \delta$ instead of $\delta$. I have a hard time understanding the latter reasoning.

The main theorem explaining the success of pre-conditioning is Theorem 3, which states that the RIP would be zero if its domain is restricted to $span_s(A)$. Theorem 4 and 5 aim to show that pre-conditioning does not deteriorate the situation for operators whose RIP was already good. So, Theorem 3 shows why pre-conditioning is expected to work and Theorem 4 and 5 show that there is no harm in using pre-conditioning. Since we never know in advance if a problem is already easy or not, we should always apply pre-conditioning and it would be detrimental if pre-conditioning makes an easy problem a hard one. We showed that this would not happen.

More precisely, we cannot claim that pre-conditioning improves the the RIP constant for nearly isometrically distributed operators. We could only claim a high probability result for the two events to happen at the same time, which are $\delta_{s}(\mathcal{A}) \leq 2\delta$ and $\delta_{s}(\mathcal{\tilde{A}}) \leq 1-(1-\delta) /[1+\sqrt{\frac{n^2}{m}}(1+\epsilon)]^2$. And as long as you can bound $\delta_{s}(\mathcal{A})$, you can also bound $\delta_{s}(\tilde {\mathcal{A}})$. Rather than saying that Gaussian (which is already an easy problem) can benefit from pre-conditioning, it is on the other side claiming that pre-conditioning will not harm even if originally the sensing operator has the largest sample efficiency.

If we look at the simulation result, we could see that for large n,m values, the empirical RIP results for original and pre-conditioned are almost the same for normal distribution and therefore there is no improvement (and no performance deterioration due to pre-conditioning). If the original sensing operator is good enough like i.i.d Gaussian, the pre-conditioning will not help much.

Dear Reviewer oF7T,

We greatly appreciate the time and effort the reviewer has put to read our work and provide constructive comments.

1. Questionable novelty

i) The idea to replace a matrix with its left singular vectors is a very natural one, and it has indeed been proposed before in the literature -- in the manuscript of Chen and Lin (2021) that the authors cite. In there, it is motivated via and formulated in the context of sparse recovery, but the mathematics are more or less equivalent. This severely impacts the value of the empirical success of their method.

The idea of preconditioning is truly not new. The main difficulty is how to transfer the low sparsity assumption from compressed sensing to the low-rank assumption in matrix sensing. In Theorems 3 and 4, we theoretically proved the RIP bound for low-rank case. By our experimental results, we show that the the preconditioning can also work in Matrix Sensing area. We also show that even for matrices with special structures like low-rank or sparse structure, pre-conditioning also works. While there are many papers on RIP, none of them has applied the pre-conditinoning idea and we believe the complexity is the analysis of low-rank matrices rather than sparse vectors. To understand the importance of this result, let us first summarize our results:

• The existing literature states that if RIP is less than 0.5, the learning problem is easy.

• Another result says that RIP can decrease below 0.5 and approach zero if the samples are from i.i.d. Gaussian and the number of samples increases.

• RIP may never become less than 1/2 for non-Gaussian even if there are infinitely many samples.

• Calculating RIP is NP-hard and it may not change continuously or smoothly with respect to changes in the probability distribution of the measurement models. Now, Section 2 asks whether RIP changes smoothly if there is a deviation from the assumptions made in the literature (i.e. from i.i.d. Gaussian).

• We showed that if the deviation from Gaussian is not too much, then increasing the number of samples would lower the RIP to the level of Gaussian and then it can become less than 0.5. To do so, we showed that there is an upper bound on RIP and it can be reduced by increasing the number of samples. This is a very strong result because for distributions far away from Gaussian it is known that even with infinitely many samples RIP can still be very large, but we showed that this is not the case if the deviation from Gaussian is modest.

• In Section 3, We consider a distribution that is far away from Gaussian.

• It is known that for such distributions, the RIP is very high and may not decrease to below 0.5 even in presence of many samples.

• We proposed a pre-conditioning technique that transforms the distribution to another one that makes the new distribution closer to Gaussian. In other words, after pre-conditioning, the individual entries of the new sensing matrices are approximately Gaussian and the preconditioned operators are likely to act as i.i.d. Gaussian with small perturbation.

• Then, the results of Section 2 would be applied to the pre-conditioned operator.

Our results allow a matrix sensting problem with many spurious solutions to be transformed to another one without any spurious solutions. This can by no means be inferred from Chen and Lin (2021) and we believe the results are powerful.

ii) I am not convinced by Lemma 1 -- the right-singular vector matrix will only be Haar distributed if the distribution of the operator A fulfills some orthogonal invariance, which is not the case for many of the distributions they use in their numerical experiments.

In "Elizabeth S. Meckes, The Random Matrix Theory of the Classical Compact Groups, Cambridge University Press, 2019" we found the following results, which support our claim that the individual entries of a random orthogonal matrix are approximately Gaussian for large matrices:

Corollary 2.6: For each $n$, let $U_n$ be a random orthogonal matrix. Then the sequence $\\{\sqrt{n}[U\_n]\_{1,1}\\}$ converges weakly to the standard Gaussian distribution, as $n \rightarrow \infty$.

Theorem 2.9: Let \left\\{U\_n\right\\} be a sequence of random orthogonal matrices with $U\_n \in \mathbb{O}(n)$ for each $n$, and suppose that $p_n q_n \xrightarrow{n \rightarrow \infty} \infty$, with $p\_n q\_n=o(n)$. Let $U\_n\left(p\_n, q\_n\right)$ denote the top-left $p_n \times q_n$ block of $U_n$, and let $Z\left(p_n, q_n\right)$ denote a $p_n \times q_n$ random matrix of i.i.d. standard normal random variables. Then

$$\lim\_{n \rightarrow \infty} d\_{TV}\left(\sqrt{n} U\_n\left(p_n, q_n\right), Z(p\_n, q\_n)\right)=0$$

Our intuition is that random orthogonal matrices have asymptotic behavior similar to Gaussian for large n values. We agree that a deeper analysis would be useful. However, all of the bounds and theoretical results we currently have in the paper are correct, and are supported by our simulations.

4. Impact of noise

• When performing the preconditioning in practice, it will not only transform the measurement operator, but also noise in measurements $y = A(x)+\epsilon$. Since the success of their method relies on the singular values of A being small (see Theorem and Remark 4), this means that the noise necessarily is amplified. A discussion on this, and possible mitigations, such as in (Chen, Lin; 2021) would make the work more complete.

Thank you for bringing up an important issue. We respectfully disagree with your conclusion that our method would amplify the noise. To understand the issue, consider the noisy case $y=\mathcal A(X)+w$, where $w$ denotes noise. The optimization problems

$

\min_X \|\mathcal A(X)-y\|\quad \text{s.t.}\quad \text{rank}(X)=r

$

and

$

\min_X \|2\mathcal A(X)-2y\|\quad \text{s.t.}\quad \text{rank}(X)=r

$

have the same solutions. This means that if we multiply all measurements by a factor of 2, the optimization solution will not change even if the noise value is increased by a factor of 2. In fact, it can be shown that the RIP value also does not change with a scalar multiplication. In fact, what matters is the relative magnitude of the elements of $\mathcal A$ compared to $w$, and multiplying both of them with the same constant will not affect their relative values. We should mention that the general definition of RIP is

$

(1-\delta)\|M\|_F^2\leq \text{constant} \times \|\mathcal A(M)\|^2\leq (1+\delta)\|M\|^2

$

where "constant" can be any arbitrary number. Many papers normalize $\mathcal A$ and consider the constant to be 1. That is why in our bounds we have $\sigma\_1(A)$, and otherwise we can re-scale it as re-scaling $\min\_X \\|\mathcal A(X)-y\\|$ will not affect the optimization solution. We will add a discussion to the paper to address your comment and emphasize that this method does not directly make the noise effect severe.

I have read the reviewers detailed response, and have carefully thought about it. Let me first say that the response 4 has helped my understanding. I was wrong to say that small singular values in of themselves are problems -- it is rather the fraction between the largest and smallest singular value (i.e. the condition number of A) that is relevant.

The rest of the responses have however not convinced me. In short, it seems to me that the authors in their response present a very different interpretation of their results in Section 3 than in the paper itself. I find this concerning for the soundness of the paper, and I will therefore not raise my score.

In their response, the authors write that

we cannot claim that pre-conditioning improves the the RIP constant for nearly isometrically distributed operators.

I.e., that the RIP-constant of the preconditioned matrix is not better than the non-preconditioned, but not significantly worse, and that that is the intention of Section 3. However, in said section ( in Remark 5), they write

Hence, the preconditioning algorithm can improve the RIP constant for nearly isometrically distributed random matrices with high probability.

These two statements diametrically contradict each other. Their updated result in particular also showcases this change of interpretation: A claim of the method improving the RIP constant to a claim that for large values of m and n, the RIP constant does not change (It should be acknowledged that their updated result does not need $m$ to grow as $n^2$, which is good).

This begs the question: Which interpretation is correct? If it is the one given in the response, I think it is fair to say that the theory does little to show that the preconditioning improves the RIP, and therefore does not explain the success of the method.

Also, although Theorem 3 gives a nice intuitive picture, it ultimately only proves that the RIP of the matrix is small if all low-rank matrices lie close to $\mathrm{span}_m(\mathcal{A})$ - whether this is true for moderately sized $m$ is left entirely open. If one relies on $\mathrm{span}_m(\mathcal{A})$ to contain the entire set of rank-$s$-matrices one will need $\mathrm{span}_m(\mathcal{A})$ to be the entire space (since $\mathrm{span}(\text{rank s matrices})=\mathbb{R}^{n,n}$), and hence that $m\geq n^2$.

Let me also clarify some things related to the responses 1.

(i) I understand that the authors formulate their results in a context of matrix sensing problems. This interpretation is nice. However, I think that as far as the theorems they prove are concerned, proving RIP for sparse signals is mathematically more or less equivalent to proving RIP for low-rank matrices (or any other set of structured signals, for that matter).

(ii) I agree that the entries of Haar matrices are approximately Gaussian distributed, I do not question the validity of this theorem. What I'm rather not convinced about is how relevant the theorem is to say something about the RIP-property of the preconditioned matrix: Unless $\mathcal{A}$ is rotationally symmetrically distributed, $\widetilde{\mathcal{A}}$ is not Haar-distributed. SInce the authors say that they want to treat distributions that are far away from Gaussian measurements (which is the example of rotationally symmetrically distributed matrices they are appeal to) in this section, it is to me not clear why Lemma 1 should apply to it.

Dear Reviewer 27fu,

We greatly appreciate the time and effort the reviewer has put to read our work and provide constructive comments.

1. the motivation of the processing algorithm is not clear: if one can design the sensing matrices at their will, what is then the advantage of starting from given ones and try to improves them? Why not taking directly Gaussian matrices with good RIP constant? This needs to be discussed at it appears nowhere.

We believe that there has been some major confusion about the main idea of our work. Before responding to your comment, we would like to provide a real-world example to better explain how $\tilde A$ is designed from $A$. The low-rank matrix sensing problem studied in this paper naturally appears in power systems, where the problem is called state estimation and is solved every 5 minutes by power system operators. A power system is a graph with $n$ nodes and a set of edges $\mathcal E$. Each node of the system has a voltage parameter $x_i$ to be learned. Each measurement $j$ of the network is in the form of

$

b_j=\sum_{i: (j,i)\in\mathcal E} \frac{x_j(x_j-x_i)}{z_{ji}}

$

where $z_{ij}$ is a known line parameter and $\frac{x_j(x_j-x_i)}{z_{ji}}$ is the power flown over line $(j,i)$. The right-hand side of the measurement $j$ can be written as

$

b_j=\langle A_j,xx^T\rangle

$

for some matrix $A_i$ that depends on the parameters $z_{ji}$ and the topology of the graph (note: $x$ is the vector of all nodal voltages). We cannot change any measurement model $A_j$ directly. Changing an entry of $A_j$ means removing/adding lines to a physical power grid or changing the reactances of the transmission lines on the streets, which is impossible (the goal is to learn the voltages from the data given by the sensors rather than changing the infrastructure). The existing methods requiring $A_j$ to be Gaussian are not applicable at all since $A_j$ is heavily structured for power systems. We propose the following idea:

• Recall that sensors $1,2,...,n$ return the measurements $b_1,...,b_n$.
• We design some coefficient $\alpha_1,...,\alpha_n$, and create a mixed measurement $\alpha_1b_1+\cdots+\alpha_nb_n$. We replace measurement 1 with this new combined measurement. Then, the new measurement can be written as $\langle \tilde A_1,xx^T\rangle$, where $\tilde A_1$ is equal to $\alpha_1 A_1+\cdots+\alpha_n A_n$.
• Note that the mixing idea cannot generate arbitrary values for $\tilde A_1$. For example, if there is no physical line between nodes 2 and 3, then the (2,3) entry of all matrices $A_1,...,A_n$ are zero and so the (2,3) entry of $\tilde A_1$ is also zero no matter what coefficients $\alpha_i$'s we select.
• We then proceed and replace measurement 2 with a new mixed measurement $\beta_1b_1+\cdots+\beta_nb_n$. We proceed with the replacement of all measurements.
• Using this idea, we exploit the existing measurements/sensors, and do not require new measurements that are not physically infeasible. The question is: how can $\alpha_i$'s, $\beta_i$'s, etc. be designed so that the process of learning $x$ becomes simpler?

To explain the above idea in a general context, the pre-conditioning trick mentioned in this paper allowed us to make linear combinations of the original sensing matrices, and we cannot change $A$ to arbitrary $\tilde{A}$ like Gaussian i.i.d. sensing operator. We use mixed measurements to obtain

$$\tilde{A}\_i = \sum\_{j=1}^m P\_{ij}A\_j,\quad \forall i\in\{1,...,m\}$$

The new sensing operator can only be in the linear span space of the original sensing matrices. And in Algorithm 1, we have

$$U,S,V^\top \leftarrow \operatorname{SVD}(\operatorname{VStack}(a_1,a_2,\ldots, a_m))$$ $$\tilde{A\_i}\leftarrow\operatorname{mat}(V\_i^\top)$$

This is equivalent to calculating a weight matrix $P = U^\top S^{-1}$ and performing the left-side multiplication $\tilde{A} = PA$. If we were to choose another matrix, say $\tilde P$, and obtain the new operator via $\tilde P A$, since $\tilde P$ has a rectangular shape, the set of possibilities for $\tilde P A$ is limited and none of the choices may be close to Gaussian. In the context of power systems, we explained that no matter how we perform the mixing, the (2,3) entry of the new operator is zero if there is no physical line (2,3) in the system.

We will revise the paper to clarify the setting of the paper in response to your constructive comment. We will include the above example to illustrate the main idea of the paper.

2. first line below def 2 is not clear. What A are we talking about that is nearly isometric? Please rephrase

Thanks. We will improve the writing. Line 118 is saying that for any A satisfying Definition 2 (Nearly isometrically distributed), the RIP constant for A will be upper bounded by $\delta$ with high probability if $m = {\Theta}(rn/\delta^2)$.

3. what is a variance "proxy" in corollary 1 or theorem 2?

Variance proxy is used in the definition of Sub-Gaussian distribution. For a sub-Gaussian random variable, there exists a parameter that bounds the tails of its distribution in a way similar to how the variance bounds the tails of a Gaussian (normal) distribution. To be more exact, if there exists some $s^2$ such that $\operatorname{E} [e^{(X-\operatorname{E}[X])t}] \leq e^{\frac{s^2t^2}{2}}$ for all $t$, then $s^2$ is called a "variance proxy".

4. rmk 2: where does it come from that A has RIP O(1/sqrt m)? Is this automatic from nearly isometrically distributed?

Yes, if $\mathcal{A}$ is nearly isometrically distributed, then there exist positive constants $c_1$ and $c_2$ such that, with probability at least $1-\exp \left(-c_1 m\right)$, we have $\delta_s(\mathcal{A}) \leq c_2 \sqrt{ns\log n /m} = \mathcal O(1/\sqrt m)$

5. the text below the proof to Theorem 3 is important but a bit too condensed to be well understood. Please expand it

Thanks. We will more exposition to the paper. For $X \in \operatorname{span}\_{s}(\mathcal A)$, we could have $\frac{||\mathcal{A}(X)||^2}{||X||\_F^2}=1$, which means that the RIP would be zero if we allow its domain to be restricted to $span\_s(A)$. As $m$ increases, the span space finally consists of all the orthonormal bases in $\mathbb{R}^{n \times n}$, and hence will cover the low-rank space $\\{X \mid \operatorname{rank}(X)\leq s\\}$ which is $\\{X\mid X=\sum\_{i=1}^s \alpha_i u_i v_i^\top\\}$. This means that the true RIP would get close to zero as we add more orthonormal sensing matrices.

For the question part:

1. In line 86, the RIP constant is defined as the smallest number such that the inequality holds; so it should be unique if it exists.

Yes, thanks for pointing out. We will change our words to "If $\delta_s$ satisfies the RIP inequality (3), every number greater than $\delta_s$ will satisfy the RIP inequality (3)".

2. In line 119, a statement is made about nearly isometry, but a distribution for A is not established. Is it the one in the following sentence? (i.i.d. Gaussian)

Not exactly; there is a wide range of distributions included in the definition of nearly isometry (examples are provided in the response to Comment 2). The statement holds for every distribution satisfying Definition 2.

3. In line 163, it appears that the mat(.) operator must know the size of the target matrix (as an additional input?)

Yes, you are correct. We need to know the size of the target matrix. Since we know the size of our sensing matrix $A$, we can do this operation. We will add this description to our updated paper.

4. In Theorem 3, is it implicitly assumed that $m>n^2$. Otherwise you would not be able to have all $A_1,\dots,A_m$ orthogonal.

No, we do not need this assumption. In fact, the problem is trivial if $m>n^2$. The practical regime is $m<n^2$ (and in fact the hope is that $m$ is close to $O(n)$ rather than $O(n^2)$ since the whole point of matrix sensing is to learn a matrix from a limited number of measurements). We could have

$$A_1 = \begin{bmatrix} 1 & 0 \\\\ 0 & 0 \\\\ \end{bmatrix}, A_2 = \begin{bmatrix} 0 & 1 \\\\ 0 & 0 \\\\ \end{bmatrix},$$

and at the same time

$$\langle A_1, A_2\rangle = \operatorname{Trace}(A_1^\top A_2)=0$$

They are orthogonal to each other, while $m =2<n^2 = 4$.

4. It is also straightforward to see that a weighted mixture of matrices will have distribution closer to Gaussian (akin to the law of large numbers). In addition, sum of subgaussian random variables are subgaussian as well (akin to Corollary 1).

There has been a major confusion here. We would like to provide a real-world example to explain our answer. The low-rank matrix sensing problem studied in this paper naturally appears in power systems, where the problem is called state estimation and is solved every 5 minutes by power system operators. A power system is a graph with $n$ nodes and a set of edges $\mathcal E$. Each node of the system has a voltage parameter $x_i$ to be learned. Each measurement $j$ of the network is in the form of

$

b_j=\sum_{i: (j,i)\in\mathcal E} \frac{x_j(x_j-x_i)}{z_{ji}}

$

where $z_{ij}$ is a known line parameter and $\frac{x_j(x_j-x_i)}{z_{ji}}$ is the power flown over line $(j,i)$. The right-hand side of the measurement $j$ can be written as

$

b_j=\langle A_j,xx^T\rangle

$

for some matrix $A_i$ that depends on the parameters $z_{ji}$ and the topology of the graph (note: $x$ is the vector of all nodal voltages). We cannot change any measurement model $A_j$ directly. Changing an entry of $A_j$ means removing/adding lines to a physical power grid or changing the reactances of the transmission lines on the streets, which is impossible (the goal is to learn the voltages from the data given by the sensors rather than changing the infrastructure). The existing methods requiring $A_j$ to be Gaussian are not applicable at all since $A_j$ is heavily structured for power systems. We propose the following idea:

• Recall that sensors $1,2,...,n$ return the measurements $b_1,...,b_n$.
• We design some coefficient $\alpha_1,...,\alpha_n$, and create a mixed measurement $\alpha_1b_1+\cdots+\alpha_nb_n$. We replace measurement 1 with this new combined measurement. Then, the new measurement can be written as $\langle \tilde A_1,xx^T\rangle$, where $\tilde A_1$ is equal to $\alpha_1 A_1+\cdots+\alpha_n A_n$.
• Note that the mixing idea cannot generate arbitrary values for $\tilde A_1$. For example, if there is no physical line between nodes 2 and 3, then the (2,3) entry of all matrices $A_1,...,A_n$ are zero and so the (2,3) entry of $\tilde A_1$ is also zero no matter what coefficients $\alpha_i$'s we select.
• We then proceed and replace measurement 2 with a new mixed measurement $\beta_1b_1+\cdots+\beta_nb_n$. We proceed with the replacement of all measurements. \item Using this idea, we exploit the existing measurements/sensors, and do not require new measurements that are not physically infeasible. The question is: how can $\alpha_i$'s, $\beta_i$'s, etc. be designed so that the process of learning $x$ becomes simpler?

As explained above, no matter how we perform the mixing, the (2,3) entry of the new operator is zero if there is no physical line (2,3) in the system. This contradicts your statement: "It is also straightforward to see that a weighted mixture of matrices will have distribution closer to Gaussian." The mixing technique does not change the number of measurements and also in the context of power systems, it cannot change the values of the zero entries corresponding to non-existent power lines, and therefore, the law of large numbers mentioned in your comment does not hold.

The underlying idea behind this confusion is similar to the issue discussed for Comment 2. $||\mathcal N||\_\infty$ being sub-Gaussian bounded does not imply i.i.d. subgaussianity of $\mathcal{N}$. We only assume that the largest magnitude of perturbation is bounded with high probability. We do not assume any distribution on perturbation, and it can be asymmetric, not centered around 0. We will add some clarification to the paper about the important issue you have raised to avoid any confusion.

5. It is not clear why Theorem 1 defines its constants in terms of $\epsilon$ when a bound of $\epsilon$ itself is already included in the assumptions; the best constant can be obtained by having $\epsilon$ meet its bound.

Thanks a lot for your advice. We will apply the upper bound of $\varepsilon$ to our result to avoid confusion. Right now the best upper bound we could have due to first-order approximation is

$$\delta\_s+ \frac{(1-\delta\_s )(1+\delta\_s)^{1/2}}{||\mathcal A||\_{\infty}}+\frac{1-\delta\_s}{8mn^2||\mathcal A||^2\_{\infty}}$$

Hence, the influence of perturbation on the RIP constant is upper bounded by a controllable constant term (We will update our paper accordingly).

2. The broadest results I am aware of for RIP of random matrices are based on subgaussianity of the underlying distribution of its entries. Thus focusing on Gaussian matrices in a comparison may be too narrow; and I note that some of the improvements claimed here involve subgaussian distributions already. Furthermore, sums of Gaussian and subgaussian matrices are subgaussian.

Thank you for your comment. The most general result of our paper is Theorem 1, which states that for arbitrary $\mathcal{A}$ and any bounded perturbation operator $\mathcal{N}$, the RIP constant will increase by at most $\mathcal{O}(mn^2\varepsilon)$. There is no assumption on subgaussianity, and this is a deterministic result.

After Theorem 1, its following Corollary and Theorem are high-probability results, based on the assumption that $\mathcal{A}$ is nearly isometrically distributed, which includes a large number of distributions. For example, it includes the i.i.d. symmetric Bernoulli distribution

$$A\_{ij}=\begin{cases} \sqrt{\frac{1}{p}} & \text { with probability } \frac{1}{2}, \\\ -\sqrt{\frac{1}{p}} & \text { with probability } \frac{1}{2},\\\ \end{cases}$$

or the following where two-thirds of the entries are zero:

$$A\_{ij}= \begin{cases}\sqrt{\frac{3}{p}} & \text { with probability } \frac{1}{6}, \\\ 0 & \text { with probability } \frac{2}{3}, \\\ -\sqrt{\frac{3}{p}} & \text { with probability } \frac{1}{6},\end{cases}$$

Regarding the deviation/perturbation, we assume $||\mathcal N||_\infty$ is sub-Gaussian bounded. Here, we need to clarify that we do not assume i.i.d subgaussianity for all entries of $\mathcal{N}$ and only require the largest magnitude of perturbation to be bounded with high probability, which is

$$\mathbf{P}(||\mathcal N||\_\infty \geq \varepsilon) \leq 2 \exp\left({-\frac{m\varepsilon^2}{\sigma^2}}\right)$$

For example, $\mathcal N$ can be [0,1] Uniform distributed with close to 1 correlation while at the same time, $||\mathcal N||\_\infty$ is upper bound by 1, hence $||\mathcal N||\_\infty$ is subgaussian. We do not assume any distribution on perturbation, and it can be asymmetric, not centered around 0, highli-correlated. Since there is no i.i.d sub-Gaussian assumption on the perturbation entries, we could not apply directly the subgaussianity result of RIP.

3. It is well known in the CS literature that orthogonalizing the rows of the sensing matrix improves its performance (these are so-called orthogonal projectors). See for example https://doi.org/10.1109/TIT.2005.862083 and https://mdav.ece.gatech.edu/publications/dwb-tr-2006.pdf - this may also relate to the Haar distributed random matrices used in Section 3.1, although they have not been defined in this paper. The preconditioning algorithm is therefore commonly applied in the CS literature, and that proposed here is a straightforward extension from measurement vectors to sensing matrices.

The idea of preconditioning is truly not new. The main difficulty is how to transfer the low sparsity assumption from compressed sensing to the low-rank assumption in matrix sensing. In Theorems 3 and 4, we theoretically proved the RIP bound for low-rank case. The existing results on orthogonalization are in a different context and cannot lead to these conclusions. We would like to emphasize that although the proposed technique may seem natural, it indeed leads to a powerful result: the RIP is often above 0.5 for which there are many spurious solution but pre-conditioning reshapes the landscape and makes the spurious solutions disappear. Many papers have been published on RIP in NeurIPS, ICML, ICLR, etc. in the past 10 years and while several of them tried to improve the RIP, none of them was able to make the observation of our paper. We believe that this is because the idea of preconditioning for low-rank matrices is not as trivial as it is for sparse vectors previously done in CS.

Dear Reviewer iS9D,

We greatly appreciate the time and effort the reviewer has put to read our work and provide constructive comments.

1. The analytical results show that the RIP constants increase with perturbation; it is not clear then what is the benefit of perturbation if it requires more measurements (e.g., Remark 2). The matrix perturbation proposed is not tested numerically.

Before responding to your comment, we would like to first summarize the goals of Section 2:

• The existing literature states that if RIP is less than 0.5, the learning problem is easy.

• Another result says that RIP can decrease below 0.5 and approach zero if the samples are from i.i.d. Gaussian and the number of samples increases.

• RIP may never become less than 1/2 for non-Gaussian even if there are infinitely many samples.

• Calculating RIP is NP-hard and it may not change continuously or smoothly with respect to changes in the probability distribution of the measurement models. Now, Section 2 asks whether RIP changes smoothly if there is a deviation from the assumptions made in the literature (i.e. from i.i.d. Gaussian).

• We showed that if the deviation from Gaussian is not too much, then increasing the number of samples would lower the RIP to the level of Gaussian and then it can become less than 0.5. To do so, we showed that there is an upper bound on RIP and it can be reduced by increasing the number of samples. This is a very strong result because for distributions far away from Gaussian it is known that even with infinitely many samples RIP can still be very large, but we showed that this is not the case if the deviation from Gaussian is modest.

Then, moving on to Section 3, we consider a distribution that is far away from Gaussian. We proposed a pre-conditioning technique that transforms the distribution to another one that makes the new distribution closer to Gaussian. Then, the results of Section 2 would be applied to the pre-conditioned operator.

In summary, we do not perturb a given operator and what we mean by perturbation is that when an original operator or a preconditioned operator is deviated from Gaussian (in which case we call it a perturbed Gaussian), we would like to understand how much the RIP would increase due to that deviation. We will improve the text of the paper and clarify that by perturbation we mean deviation.

In Section 2, we derive the upper bound for RIP under deviation from (perturbation to) Gaussian. In Theorem 2, based on our assumption, the perturbation is mean 0 and variance proxy $\sigma^2/m$ (to match with the magnitude "Gaussian i.i.d. sensing operators $N(0,1/m)$ satisfy RIP"), and that for $\sigma$

$$\sigma<\frac{1-\delta\_s}{2c\sqrt{m}n^2\|\mathcal A\|\_{\infty}}$$

The RIP upper bound we derive is $\mathcal{O}(1/\sqrt{m})+\mathcal{O}(\sqrt{m}\sigma)+\mathcal{O}(\sigma^2)$, and hence the influence of perturbation on RIP is a controllable constant level. This also means that we can compensate for the RIP increase due to deviation from Gaussian by slightly increasing the number of measurements $m$. And that is the intuition behind the pre-conditioning algorithm. After pre-conditioning, the individual entries of the new sensing matrix are approximately Gaussian, and hence these preconditioned operators are likely to act as i.i.d. Gaussian with small perturbation. As long as the new upper bound after perturbation is smaller than 0.5, we can obtain good properties like global optimality for Matrix Sensing problems. We will improve the presentation of our results based on your constructive comment.

For the numerical result of matrix perturbation, we could see that there is a small influence on the empirical RIP value as long as the assumption for perturbation holds: $\sigma<\frac{1-\delta_s}{2c\sqrt{m}n^2\|\mathcal A\|_{\infty}}$.

Thanks for the response. My interpretation of N as a perturbation of A was based on the context given by the paper title (“manipulation”) and the perturbation itself (which is given first as an addition of a deterministic bounded matrix, instead of described in terms of the distribution of A). I can conceive of a practical setting for this setup, e.g., when the matrix is generated according to a Gaussian distribution and is then quantized. The second setting is more difficult to assess in terms of practical impact, as both A and N are drawn from specific distributions and, as other reviewers have pointed out, it’s not clear why it is better to perturb and add measurements than to keep the original matrix distribution at the larger size, which should also improve the RIP constant.

Also, my last question was wrong and the matrices can only be orthogonal if m < n^2, so we are in agreement.

2. The justification of the pre-conditioning is not convincing enough. In Section 2, they prove the RIP constant will increase in the presence of perturbation. There are several issues. First, their analysis is on the upper-bound on RIP. It's unclear how the RIP constant change. The RIP constant can remain constant while the upper bound of RIP constant can change from c to 2c. Second, the upper-bound is a little confusing for me. It seems that there a optimal sensing number m. Intuitively, the performance should keep improving with the increasing sampling number, as the deviation becomes smaller.

Section 2 is regarding perturbation analysis, while Section 3 is about pre-conditioning. Let us first summarize the goals of Section 2:

• The existing literature states that if RIP is less than 0.5, the learning problem is easy.

• Another result says that RIP can decrease below 0.5 and approach zero if the samples are from i.i.d. Gaussian and the number of samples increases.

• RIP may never become less than 1/2 for non-Gaussian even if there are infinitely many samples. So, your statement that the perform should improve (or RIP should decrease) by having more samples is unfortunately proven not to hold in general.

• Calculating RIP is NP-hard and it may not change continuously or smoothly with respect to changes in the probability distribution of the measurement models. Now, Section 2 asks whether RIP changes smoothly if there is a deviation from the assumptions made in the literature (i.e. from i.i.d. Gaussian).

• We showed that if the deviation from Gaussian is not too much, then increasing the number of samples would lower the RIP to the level of Gaussian and then it can become less than 0.5. To do so, we showed that there is an upper bound on RIP and it can be reduced by increasing the number of samples. This is a very strong result because for distributions far away from Gaussian it is known that even with infinitely many samples RIP can still be very large, but we showed that this is not the case if the deviation from Gaussian is modest.

Now, let us summarize the goals of Section 3:

• We consider a distribution that is far away from Gaussian.
• It is known that for such distributions, the RIP is very high and may not decrease to below 0.5 even in presence of many samples.
• We proposed a pre-conditioning technique that transforms the distribution to another one that makes the new distribution closer to Gaussian. In other words, after pre-conditioning, the individual entries of the new sensing matrices are approximately Gaussian and the preconditioned operators are likely to act as i.i.d. Gaussian with small perturbation.
• Then, the results of Section 2 would be applied to the pre-conditioned operator.

Regarding your second concern, that is due to technical details of the proof. We also assume that the perturbation is mean 0 and variance proxy $\sigma^2/m$, so as to match with the magnitude "Gaussian i.i.d. sensing operators $N(0,1/m)$ satisfy RIP". Our RIP bound in Theorem 2 can be decomposed into two parts: the first is the RIP bound for $\mathcal{A}$ to be nearly isometrically distributed, and the second part is influence by noise. As m increases, the first part decreases, and the second part seemingly increases on the order of $\mathcal{O}(\sqrt{m}\sigma)$. However, we could see that there is an assumption for $\sigma$:

$$\sigma<\frac{1-\delta\_s}{2c\sqrt{m}n^2\|\mathcal A\|\_{\infty}}$$

By taking the magnitude of $\sigma$ into consideration, your intuition is absolutely correct. The performance keeps improving with increasing the number of samples, as the deviation becomes smaller.

3. The techniques in the analysis is quite standard.

The problem defined in this paper is new to Matrix Sensing, which is how to improve the landscape when it is not possible to change the sensors that collect data (as discussed through the power systems example above). The problem has previously been studied in signal processing for sparse vectors. However, in our paper, the analysis is for low-rank matrices. The framework is similar but details differ. Although the technique may seem standard, the observation that a simple mixing technique makes the RIP become less than 0.5 is powerful and very useful for practical purposes.

4. A minor regarding to Figure 2 is to plot the error bar of the simulated RIP constant. Unless you exhaustively check all possible rank-s matrices, there is no possibility of getting the exact RIP constant. Based on the content, it seems that you use Monte-Carlo simulation to find the RIP constant, which is fine. Still, it looks better to include the variance.

Thanks a lot for your advise. We have re-run the experiments and add variance to the plot.

Dear Reviewer aJfp,

We greatly appreciate the time and effort the reviewer has put to read our work and provide constructive comments.

1. My major concern comes from the setting itself. The authors basically change the sensing operator, e.g., from $A$ to $\tilde{A}$ in Algorithm 1. If this operation is allowed, I wonder why the authors bother with the pre-conditioning trick, they can simply ignore the original operator A and use a Gaussian i.i.d. sensing operator. The RIP constant can be guaranteed.

We believe that there has been some major confusion about the main idea of our work. Before responding to your comment, we would like to provide a real-world example to better explain how $\tilde A$ is designed from $A$. The low-rank matrix sensing problem studied in this paper naturally appears in power systems, where the problem is called state estimation and is solved every 5 minutes by power system operators. A power system is a graph with $n$ nodes and a set of edges $\mathcal E$. Each node of the system has a voltage parameter $x_i$ to be learned. Each measurement $j$ of the network is in the form of

$

b_j=\sum_{i: (j,i)\in\mathcal E} \frac{x_j(x_j-x_i)}{z_{ji}}

$

where $z_{ij}$ is a known line parameter and $\frac{x_j(x_j-x_i)}{z_{ji}}$ is the power flown over line $(j,i)$. The right-hand side of the measurement $j$ can be written as

$

b_j=\langle A_j,xx^T\rangle

$

for some matrix $A_i$ that depends on the parameters $z_{ji}$ and the topology of the graph (note: $x$ is the vector of all nodal voltages). We cannot change any measurement model $A_j$ directly. Changing an entry of $A_j$ means removing/adding lines to a physical power grid or changing the reactances of the transmission lines on the streets, which is impossible (the goal is to learn the voltages from the data given by the sensors rather than changing the infrastructure). The existing methods requiring $A_j$ to be Gaussian are not applicable at all since $A_j$ is heavily structured for power systems. We propose the following idea:

• Recall that sensors $1,2,...,n$ return the measurements $b_1,...,b_n$.
• We design some coefficient $\alpha_1,...,\alpha_n$, and create a mixed measurement $\alpha_1b_1+\cdots+\alpha_nb_n$. We replace measurement 1 with this new combined measurement. Then, the new measurement can be written as $\langle \tilde A_1,xx^T\rangle$, where $\tilde A_1$ is equal to $\alpha_1 A_1+\cdots+\alpha_n A_n$.
• Note that the mixing idea cannot generate arbitrary values for $\tilde A_1$. For example, if there is no physical line between nodes 2 and 3, then the (2,3) entry of all matrices $A_1,...,A_n$ are zero and so the (2,3) entry of $\tilde A_1$ is also zero no matter what coefficients $\alpha_i$'s we select.
• We then proceed and replace measurement 2 with a new mixed measurement $\beta_1b_1+\cdots+\beta_nb_n$. We proceed with the replacement of all measurements.
• Using this idea, we exploit the existing measurements/sensors, and do not require new measurements that are not physically infeasible. The question is: how can $\alpha_i$'s, $\beta_i$'s, etc. be designed so that the process of learning $x$ becomes simpler?

To explain the above idea in a general context, the pre-conditioning trick mentioned in this paper allowed us to make linear combinations of the original sensing matrices, and we cannot change $A$ to arbitrary $\tilde{A}$ like Gaussian i.i.d. sensing operator. We use mixed measurements to obtain

$$\tilde{A}_i =\sum\_{j=1}^{m} P\_{ij} A_j, \quad \forall i \in \{1,...,m\}$$

The new sensing operator can only be in the linear span space of the original sensing matrices. And in Algorithm 1, we have

$$U,S,V^\top \leftarrow \operatorname{SVD}(\operatorname{VStack}(a_1,a_2,\ldots, a_m))$$ $$\tilde{A_i}\leftarrow\operatorname{mat}(V_i^\top)$$

This is equivalent to calculating a weight matrix $P = U^\top S^{-1}$ and performing the left-side multiplication $\tilde{A} = PA$. If we were to choose another matrix, say $\tilde P$, and obtain the new operator via $\tilde P A$, since $\tilde P$ has a rectangular shape, the set of possibilities for $\tilde P A$ is limited and none of the choices may be close to Gaussian. In the context of power systems, we explained that no matter how we perform the mixing, the (2,3) entry of the new operator is zero if there is no physical line (2,3) in the system.

We will revise the paper to clarify the setting of the paper in response to your constructive comment. We will include the above example to illustrate the main idea of the paper.