Low-Rank Matrix Approximation in the Infinity Norm

TL;DR

This paper investigates the computational complexity of low-rank matrix approximation under the entry-wise infinity norm, proving NP-completeness for rank-1 cases, analyzing solvable cases, and proposing a heuristic for practical applications.

Contribution

It establishes NP-completeness for rank-1 approximation, analyzes polynomial-time solvable cases, and introduces a heuristic algorithm for quantized matrix recovery.

Findings

Rank-1 approximation decision problem is NP-complete.

Certain cases of the problem are solvable in polynomial time.

A heuristic algorithm effectively recovers quantized low-rank matrices.

Abstract

The low-rank matrix approximation problem with respect to the entry-wise $\ell_{\infty}$-norm is the following: given a matrix $M$ and a factorization rank $r$, find a matrix $X$ whose rank is at most $r$ and that minimizes $\max_{i,j} |M_{ij} - X_{ij}|$. In this paper, we prove that the decision variant of this problem for $r=1$ is NP-complete using a reduction from the problem `not all equal 3SAT'. We also analyze several cases when the problem can be solved in polynomial time, and propose a simple practical heuristic algorithm which we apply on the problem of the recovery of a quantized low-rank matrix.

Tables

Table 1. Table 1 : Repartition of the number of iterations performed by the secant method to solve ( 4.1 ) among 10 4 superscript 10 4 10^{4} instances and for different values of m 𝑚 m , generating each entry of M 𝑀 M and u 𝑢 u using the normal distribution N ​ ( 0 , 1 ) 𝑁 0 1 N(0,1) .

m / # it.	1	2	3	4	5	6	7	8	9	Time (s.)
10	2019	1594	3277	2742	363	5	0	0	0	5.1
${10}^2$	199	204	487	4237	4099	747	27	0	0	6.55
${10}^3$	20	19	50	2668	5200	1885	154	4	0	8.13
${10}^4$	2	1	11	1739	5118	2769	351	9	0	22.03
${10}^5$	0	0	0	1258	4785	3403	525	28	1	108.11
${10}^6$	0	0	0	916	4618	3714	712	38	2	1685.40
${10}^7$	0	0	0	687	4192	4115	945	60	1	16793.87

Table 2. Table 2: Results of Algorithm  1 on Gaussian random instances of the recovery of quantized low-rank matrices.

$r$	Error of Alg. 1	# runs	# it. of Alg. 1	Average	Error of ${\ell }_2$ LRA (init.)
	min , mean , max	error $\le 0.5$	min , mean , max	time (s.)	min , mean , max
1	0.50 , 0.50 , 0.50	100/100	14 ,  49 , 110	0.67	0.92 , 0.96 , 0.97
2	0.53 , 0.69 , 0.87	0/100	3 ,   5 ,  36	0.16	0.82 , 0.94 , 0.97
5	0.52 , 0.54 , 0.62	0/100	3 ,  19 , 127	1.47	0.65 , 0.70 , 0.87
10	0.50 , 0.52 , 0.56	0/100	5 ,  90 , 265	15.19	0.70 , 0.75 , 0.85
20	0.48 , 0.49 , 0.53	93/100	14 , 175 , 333	58.69	0.74 , 0.81 , 0.91

Table 3. Table 3: Results for the BCD Algorithm  1 initialized using the solution to ℓ 2 subscript ℓ 2 \ell_{2} LRA of matrix M q subscript 𝑀 𝑞 M_{q} (exactly as for the second and last row of Table  2 ) and with the solution to ℓ 2 subscript ℓ 2 \ell_{2} LRA of M = U ​ V 𝑀 𝑈 𝑉 M=UV , where each entry of U 𝑈 U and V 𝑉 V is generated using the Gaussian distribution N ​ ( 0 , 1 ) 𝑁 0 1 N(0,1) .

	Error of BCD	# runs	# it. of BCD	Average
	min , mean , max	error $\le 0.5$	min , mean , max	time (s.)
$r=2$, init. ${\ell }_2$	0.52 , 0.67 , 0.88	0/100	3 , 6 , 35	0.19
$r=2$, init. $M$	0.50 , 0.50 , 0.50	100/100	2 , 10 , 51	0.31
$r=20$, init. ${\ell }_2$	0.48 , 0.49 , 0.51	96/100	71 , 178 , 261	56.12
$r=20$, init. $M$	0.46 , 0.46 , 0.46	100/100	66 , 96 , 135	30.60

Equations

$$\underset{X\in \Omega }{\min }||M-X||\text{ such that }\mathrm{rank}(X)\le r,$$

$$\underset{X\in {\mathbb{R}}^{m\times n}}{\min }||M-X|{|}_{\infty }\text{ such that }\mathrm{rank}(X)\le r,$$

$$M=\left(\begin{array}[]{cc}1&1\\ 1&-1\end{array}\right).$$

$$\underset{\mathrm{rank}(X)\le 1}{\min }||M-X|{|}_{\infty }=\underset{u\in {\mathbb{R}}^2,v\in {\mathbb{R}}^2}{\min }||M-u{v}^T|{|}_{\infty }=||M|{|}_{\infty }=1.$$

$$-k\le |M_{ij}-u_i{v}_j|\le k⟺M_{ij}-k\le u_i{v}_j\le M_{ij}+k.$$

$$s_i(M_{ij}-k)\le v_j\le s_i(M_{ij}+k)\text{ for all }i,j.$$

$$|M_{ij}|>k\Rightarrow \mathrm{sign}(u_i{v}_j)=\mathrm{sign}(M_{ij}),$$

$$\operatorname{sign}(x)=\left\{\begin{array}[]{cc}1&\text{ if $x>0$,}\\ 0&\text{ if $x=0$,}\\ -1&\text{ if $x<0$.}\end{array}\right.$$

$$u_i=\frac{1}{\sqrt{2}}-i\epsilon ,v_j=\frac{1}{\sqrt{2}}+j\epsilon +{\epsilon }^{1.5},$$

$$M=\left(\begin{array}[]{ccccc}2&0&1&1&-1\\ -1&2&-1&-1&0\\ -1&1&2&-1&-1\\ -1&1&1&2&-1\\ 1&-1&0&1&2\\ \end{array}\right),$$

$$N=\left(\begin{array}[]{ccccc}2&1&1&1&0\\ 0&2&1&1&1\\ -1&0&2&1&1\\ -1&-1&-1&2&1\\ -1&-1&-1&-1&2\\ \end{array}\right).$$

$$M=\left(\begin{array}[]{ccccc}2&1&1&-1&1\\ -1&2&-1&-1&0\\ -1&0&2&1&1\\ 0&1&-1&2&-1\\ -1&-1&-1&1&2\\ \end{array}\right)$$

$$\underset{v_j}{\min }\underset{i}{\max }|M_{ij}-u_i{v}_j|.$$

$$\underset{v_j}{\min }\underset{\{i|u_i\ne 0\}}{\max }|M_{ij}-u_i{v}_j|,$$

$$v_{j}^{*}\in\left\{\frac{M_{i_{1}j}+M_{i_{2}j}}{u_{i_{1}}+u_{i_{2}}}\ \big{|}\ 1\leq i_{1}\neq i_{2}\leq m,u_{i_{1}}>0,u_{i_{2}}>0\right\}.$$

$$i_l={\mathrm{argmin}}_{\{i|u_i\ne 0\}}\frac{M_{ij}}{u_i}\text{ and }{i}_u={\mathrm{argmax}}_{\{i|u_i\ne 0\}}\frac{M_{ij}}{u_i}.$$

$$M=\left(\begin{array}[]{cccccccc}0.35&0.65&0.15&0.54&1.49\\ 1.17&-0.90&-1.50&-0.52&-0.44\\ 1.03&-1.12&-3.48&-1.41&-0.17\\ 4.46&-1.81&3.58&2.24&-2.23\\ -1.53&-0.60&-2.73&-1.94&-1.10\\ -2.53&2.79&1.02&0.75&3.59\\ 3.38&-0.90&-1.16&0.53&0.86\\ -0.66&0.42&0.71&0.20&0.14\\ \end{array}\right)$$

$$M_{q}=\left(\begin{array}[]{cccccccc}0&1&0&1&1\\ 1&-1&-1&-1&0\\ 1&-1&-3&-1&0\\ 4&-2&4&2&-2\\ -2&-1&-3&-2&-1\\ -3&3&1&1&4\\ 3&-1&-1&1&1\\ -1&0&1&0&0\\ \end{array}\right),$$

$$\{X|\mathrm{rank}(X)=3,||M-X|{|}_{\infty }\le 1/2\}$$

Full text

Keywords. low-rank matrix approximations, $\ell_{\infty}$ norm, computational complexity

1 Introduction

Low-rank matrix approximations (LRA) are key problems in numerical linear algebra, and have become a central tool in data analysis and machine learning; see, e.g., [21]. A possible formulation for LRA is the following: given a matrix $M\in\mathbb{R}^{m\times n}$ and a factorization rank $r$, solve

[TABLE]

for some given (pseudo) norm $||.||$. In this paper, we focus on unconstrained variants, that is, $\Omega=\mathbb{R}^{m\times n}$, although there exists many important variants of (1.1) that take constraints into account, e.g., nonnegative matrix factorization [16], independent component analysis [4] and sparse principal component analysis [5].

When the norm $||.||$ is the Frobenius norm, that is, $||M-X||_{F}^{2}=\sum_{i,j}(M_{ij}-X_{ij})^{2}$, the problem can be solved using the singular value decomposition and is closely related to principal component analysis [10]. In practice, it is often required to use other norms, e.g., the $\ell_{1}$ norm which is more robust to outliers [20], weighted norms that can be used when data is missing [6, 15], $\sum_{j}||M(:,j)-X(:,j)||_{2}^{p}$ for $p\geq 1$ which can model different situations depending on the value of $p$ [3], and Kullback-Leibler (KL) divergence when Poisson noise is present [1]. However, as soon as the norm is not the Frobenius norm, (1.1) becomes difficult in general; in particular it was proved to be NP-hard for all the previously listed cases except for the KL divergence [9, 8, 20]. (For the KL divergence, proving NP-hardness is, to the best of our knowledge, an open problem.)

In this paper, we focus on the variant with the component-wise $\ell_{\infty}$ norm:

[TABLE]

where $||M-X||_{\infty}=\max_{i,j}|M_{ij}-X_{ij}|$. We will refer to this problem as $\ell_{\infty}$ LRA. It should be used when the noise added to the low-rank matrix follows an i.i.d. uniform distribution. We will also use the notation $\ell_{p}$ LRA for the LRA problem where the norm used is the component-wise $\ell_{p}$ norm.

1.1 Previous results and applications

When $m=n$ and $r=\min(m,n)-1$, (1.2) was studied in [18] and corresponds to the problem of distance to robust nonsingularity which can be stated as follows: what is the rank-deficient matrix $X$ that is the closest to $M$ in the component-wise infinity norm? In this particular case, (1.2) was shown to be NP-hard [18]. In [11, 12], Goreinov and Tyrtyshnikov obtained an approximate solution to $\ell_{\infty}$ LRA using the so-called rank-$r$ skeleton approximation which uses $r$ columns and $r$ rows of the given matrix. Juditsky et al. [14] linked a variant of (1.2) (where the row range of $X$ is constrained to be contained in the row range of $M$) with the synthesis problem of compressed sensing and provided randomized algorithms with optimality guarantees. Very recently, Chierichetti et al. [2] proposed provably good approximation algorithms for $\ell_{p}$ LRA for any $p\geq 1$ using a subset of the columns of the input matrix $M$ to span the columns of $X$. An application of (1.2) is the recovery of a low-rank matrix from a quantization [13] (roughly speaking, its rounding to some precision); see, e.g., [7] for more details on quantization. For example, assume we are given a real rank-$r$ matrix where each entry has been rounded to the nearest integer. Given such a matrix $M$ (which in general will have full rank), the problem is to find a rank-$r$ $X$ such that $||M-X||_{\infty}\leq 0.5$ (note that rounding can be assimilated to noise distributed uniformly in the interval [-0.5,0.5]); see Section 4 for some examples.

As far as we know, there is not as much literature on $\ell_{\infty}$ LRA (1.2) compared to other variants. A plausible explanation is that, in practice, and especially in data analysis applications, using the $\ell_{\infty}$ norm is not very useful and is in particular extremely sensitive to outliers. Moreover, even if $M\neq 0$, the zero matrix could be an optimal solution of (1.2) (which is not possible for any other $\ell_{p}$-norm).

Example 1**.**

Let

[TABLE]

We have

[TABLE]

In fact, assume the minimum is strictly smaller than 1. This requires $u_{1}v_{i}>0$ $(i=1,2$), $u_{2}v_{1}>0$ and $u_{2}v_{2}<0$ which is not possible.

However, there are several applications as mentioned above; namely, distance to singularity, compressed sensing, and recovery of a quantized low-rank matrix. Except for the case $m=n$ and $r=n-1$ which was proved to be NP-hard by Poljak and Rohn [18], there is, to the best of our knowledge, not a very good understanding of the computational complexity of (1.2). The main goal of this paper is to shed some light on this question; in particular proving NP-completeness of (1.2) when $r=1$ (Theorem 3).

1.2 Outline of the paper and contributions

In this paper we mainly focus on rank-one $\ell_{\infty}$ LRA, that is, (1.2) with $r=1$. In Section 2, we show that the decision version of rank-one $\ell_{\infty}$ LRA is in NP. In fact, we show that if the sign pattern of $X$ is known, then the decision version of rank-one $\ell_{\infty}$ LRA can be solved in polynomial time by finding a solution to a system of linear inequalities. This is an important result because it turns out that, in most cases, the number of possible sign patterns of $X$ can be reduced drastically; e.g., $X$ can be assumed to be nonnegative if $M$ is. In Section 3, we prove that rank-one $\ell_{\infty}$ LRA is NP-complete using a reduction from ‘not all equal 3SAT’, with an intermediate problem on directed graphs (namely, the problem of making, if possible, a directed graph acyclic by reversing the direction of a particular subset of the edges). Finally, in Section 4, we propose a simple heuristic algorithm and apply it on the recovery of quantized low-rank matrices.

2 Decision version of rank-one $\ell_{\infty}$ LRA

Let us define formally the decision version of rank-one $\ell_{\infty}$ LRA:

Problem 1**.**

(D-$\ell_{\infty}$-R1A($M$,$k$))

Given: A real $m$-by-$n$* matrix $M$ and a real number $k\geq 0$.*

Question: Does there exist $u\in\mathbb{R}^{m}$ and $v\in\mathbb{R}^{n}$ such that $\max_{i,j}|M_{ij}-u_{i}v_{j}|\leq k$? If yes, output a solution $(u,v)$.

Lemma 1**.**

If the sign pattern of $u$ or $v$ is given as a part of the data, then D-$\ell_{\infty}$-R1A($M$,$k$) can be solved in polynomial time in $m$ and $n$, namely in O\big{(}mn(m+n)^{3}\log mn\big{)} arithmetic operations.

Proof.

First, we can assume without loss of generality (w.l.o.g.) that for each row (resp. column) of $M$, there is at least one entry whose absolute value is larger than $k$, that is, for all $i$ (resp. $j$), there exists $l$ (resp. $p$) such that $|M_{il}|>k$ (resp. $|M_{pj}|>k$). In fact, if it is not the case, then one can trivially choose $u_{i}=0$ (resp. $v_{j}=0$) in a solution of D-$\ell_{\infty}$-R1A($M$,$k$) and reduce the problem to a submatrix of $M$. This implies that we can assume w.l.o.g. that $u\neq 0$ and $v\neq 0$ since $|M_{ij}|>k\Rightarrow u_{i}v_{j}\neq 0$.

Second, if we assume that the sign pattern of $u$ is known we can assume w.l.o.g. that $u>0$. In fact, if some entries of $u$ are negative, we can flip their signs along with the signs of the entries of the corresponding rows of $M$ to obtain an equivalent problem.

Therefore, if the sign pattern of $u$ is known, D-$\ell_{\infty}$-R1A($M$,$k$) can be reduced to finding $u>0$ and $v$ such that

[TABLE]

Moreover, by the scaling degree of freedom (that is, $uv^{T}=(\alpha u)(\alpha^{-1}v)^{T}$ for any $\alpha\neq 0$), we can assume w.l.o.g. that $u\leq 1$. Defining $s_{i}=(u_{i})^{-1}\geq 1$, this problem is equivalent to finding $s\geq 1$ and $v$ such that

[TABLE]

This is a system of $2mn+m$ linear inequalities with $m+n$ variables, where each inequality contains at most two variables, and can be solved in O\big{(}IN^{3}\log I\big{)}=O\big{(}mn(m+n)^{3}\log mn\big{)} arithmetic operations where $I=2mn+m$ is the number of inequalities and $N=m+n$ the number of variables [17].

If the sign pattern of $v$ is known, by symmetry ($||M-uv^{T}||_{\infty}=||M^{T}-vu^{T}||_{\infty}$), the same result holds. This completes the proof. ∎

Corollary 1**.**

If $M\geq 0$, D-$\ell_{\infty}$-R1A($M$,$k$) can be solved in polynomial time.

Proof.

In fact, in that case, one can assume w.l.o.g. that $u\geq 0$ and $v\geq 0$ so that the result follows from Lemma 1. ∎

Remark 1**.**

It is interesting to note that some rank-one LRA problems are NP-hard even when $M\geq 0$; e.g., for the $\ell_{1}$ norm [9] and weighted norms [8].

Lemma 1 implies that D-$\ell_{\infty}$-R1A($M$,$k$) can be solved in O\big{(}2^{\min(m,n)}mn(m+n)^{3}\log mn\big{)} operations since one can try all the possible sign patterns for $u$ (or $v$ if $n\leq m$). It is possible to achieve a better complexity result by identifying the connected component of a particular bipartite graph.

Definition 1**.**

Given $M$ and $k$, $G_{b}(M,k)=(V_{1}\cup V_{2},E)$ is the bipartite graph with $V_{1}=\{v_{1},v_{2},\dots,v_{m}\}$, $V_{2}=\{v_{1}^{\prime},v_{2}^{\prime},\dots,v_{n}^{\prime}\}$, and $(v_{i},v_{j}^{\prime})\in E\iff|M_{ij}|>k$.

Lemma 2**.**

For D-$\ell_{\infty}$-R1A($M$,$k$), the number of possible sign patterns in a solution $u$ can be reduced to $2^{d-1}$ where $d$ is the number of connected components of $G_{b}(M,k)$ discarding the isolated vertices.

Proof.

Note that the isolated vertices of $G_{b}(M,k)$ correspond to rows and columns of $M$ whose entries have absolute value smaller than $k$ and for which one can set w.l.o.g. the corresponding entry of $u$ or $v$ to zero; see the proof of Lemma 1.

The rest of the proof follows from the fact that any solution $(u,v)$ of D-$\ell_{\infty}$-R1A($M$,$k$) must satisfy the following:

[TABLE]

where

[TABLE]

Therefore, fixing the sign of an entry of $u$ imposes the sign for all the entries of $u$ and $v$ contained in the same connected component. This makes $2^{d}$ possible sign patterns for $u$ and $v$. It can be reduced to $2^{d-1}$ since $uv^{T}=(-u)(-v)^{T}$ hence half of the possible sign patterns can be discarded (this can be achieved for example by imposing arbitrarily the sign of one entry of $u$ or $v$). ∎

For any submatrix of $M$ corresponding to a connected component of $G_{b}(M,k)$, there must exist a completion with $\pm 1$ of the entries smaller than $k$ such that its sign pattern has rank one, otherwise the answer to D-$\ell_{\infty}$-R1A($M$,$k$) is NO; see, e.g., Example 1 for any $k<1$. This observation can be used to quickly obtain a lower bound on $k$ in order for the answer to D-$\ell_{\infty}$-R1A($M$,$k$) to possibly be YES (hence also a lower bound for rank-one $\ell_{\infty}$ LRA).

Theorem 1**.**

D-$\ell_{\infty}$-R1A($M$,$k$) is in NP, and can be solved in O\big{(}2^{d}mn(m+n)^{3}\log mn\big{)} operations where $d$ is the number of connected components of $G_{b}(M,k)$ discarding isolated vertices.

Proof.

This follows directly from Lemmas 1 and 2. ∎

Theorem 1 implies that for D-$\ell_{\infty}$-LRA to be a difficult problem, the number of connected components has to be high. This will motivate our construction in our NP-completeness proof where we will use a square matrix for which only the diagonal entries are larger than $k$ so that the number of connected components is maximal, namely $d=m=n$.

Theorem 1 also implies that D-$\ell_{\infty}$-R1A($M$,$k$) can be solved in polynomial time if the number of connected components $d$ satisfies d=O\big{(}\log(\min(m,n))\big{)}.

3 NP-completeness of D-$\ell_{\infty}$-R1A($M$,$k$)

The goal of this section is to prove that D-$\ell_{\infty}$-R1A is NP-hard. In order to do this, we construct a polynomial time reduction from the problem known as NOT-ALL-EQUAL 3-SAT. Recall that a literal associated with a set $X$ of Boolean variables is either an element of $X$ or a negation of it.

Problem 2**.**

(NOT-ALL-EQUAL $3$-SAT)

Given: A set $X$ of variables and a set $L$ of $3$-tuples of literals.

Question: Does there exist an assignment of the variables in $X$ to $\{0,1\}$ for which every tuple in $L$ has at least one false literal and at least one true literal?

Since NOT-ALL-EQUAL $3$-SAT is NP-complete (see [19]), constructing a polynomial time reduction from it to D-$ell_{\infty}$-R1A would mean the NP-hardness of the latter problem. In order to present such a reduction, we need to recall the definitions of some basic concepts in graph theory. An oriented graph $G$ is defined as a pair of sets $V$ and $E\subset V^{2}$. The elements of $V$ are called vertices, a pair $(a,b)\in E$ is an edge passing from $a$ to $b$, and vertices $a,b$ are adjacent if there is an edge passing between them. We assume that $V=\{1,\ldots,n\}$ and that at most one of the pairs $(a,b)$, $(b,a)$ belongs to $E$ for all $a,b\in V$. A sequence $(a_{0},\ldots,a_{k})$ of vertices is called a cycle if $a_{0}=a_{k}$ and there is an edge passing from $a_{i-1}$ to $a_{i}$ for all $i$. A two-coloring of $G$ is a partition of $V$ into the union of two disjoint sets $W$ and $B$. A subset $U\subset V$ is called monochromatic with respect to $(W,B)$ if either $U\subset W$ or $U\subset B$. Let us introduce the auxiliary problem which we use as a tool in our NP-hardness proof.

Problem 3**.**

Given: An oriented graph $G=(V,E)$ and a set $D$ of pairs of non-adjacent vertices.

Question: Is there a two-coloring $(W,B)$ of $V$ such that

• (i)

no pair in $D$ is monochromatic, and

• (ii)

the graph obtained from $G$ by reversing the edges passing between $W$ and $B$ has no cycle.

Lemma 3**.**

Problem 3 is NP-complete.

Proof.

Let us construct the graph $G$ depending on an instance $(X,L)$ of Problem 2:

• Step 1.

For every variable $x\in X$, we create two vertices corresponding to $x$ and the negation of $x$, and we add to $D$ the pair containing these two.

• Step 2.

For every tuple $(y_{1},y_{2},y_{3})\in L$, we create three new vertices corresponding to $y_{1}$, $y_{2}$ and $y_{3}$, and we draw a cycle on these vertices.

• Step 3.

We add a pair in $D$ containing a vertex created in Step 1 with a vertex created in Step 2 if they correspond to literals that are negations of each other.

Clearly, this graph can be constructed in polynomial time, and a two-coloring of $V$ leaves no pair in $D$ monochromatic if and only if it assigns the same color for every occurrence of a literal $y$ and the other color for the negation of $y$. If $y_{1}$, $y_{2}$ and $y_{3}$ have all the same color and $(y_{1},y_{2},y_{3})\in L$, then the item (ii) in Problem 3 does not require changes in the edge directions between $y_{1}$, $y_{2}$ and $y_{3}$, so these vertices remain a cycle. This implies that any acceptable two-coloring of $G$ for Problem 3 will not have $y_{1}$, $y_{2}$ and $y_{3}$ of the same color hence will correspond to a valid assignment for $(X,L)$. On the other hand, any valid assignment of $(X,L)$ corresponds to a coloring in which $y_{1}$, $y_{2}$ and $y_{3}$ have different colors for all $(y_{1},y_{2},y_{3})\in L$. In other words, two of the three edges passing between vertices in $y_{1}$, $y_{2}$ and $y_{3}$ will change their directions as in item (ii) in Problem 3, which means that the resulting graph will possess no cycle hence any valid assignment of $(X,L)$ corresponds to an acceptable two-coloring of $G$. ∎

We are now ready to present a reduction from Problem 3 to D-$ell_{\infty}$-R1A.

Definition 2**.**

Let $G=(V,E)$ and $D$ be defined as in the formulation of Problem 3. We define the matrix $\mathcal{M}=\mathcal{M}(G,D)\in\{-1,0,1\}^{|V|\times|V|}$ with rows and columns indexed by elements of $V=\{1,2,\dots,n\}$ as follows:

• (1)

$\mathcal{M}_{ii}=2$* for all $i$,*

• (2)

$\mathcal{M}_{ij}=\mathcal{M}_{ji}=-1$* if $\{i,j\}\in D$,*

• (3)

$\mathcal{M}_{ij}=-1$, $\mathcal{M}_{ji}=1$ if $(i,j)\in E$,

• (4)

$\mathcal{M}_{ij}=0$* otherwise.*

Before we can prove that $\mathcal{M}$ leads to a reduction, we need a result stating that any partial order relation on a set can be extended to a total order relation. In terms of graphs, this result can be stated as follows.

Observation 1**.**

Let $G=(V,E)$ be an oriented graph without cycles. Then there exists a total ordering $\succ$ of $V$ such that $(a,b)\in E$ implies $a\succ b$.

Theorem 2**.**

The pair $(G,D)$ is a yes-instance of Problem 3 if and only of there are real numbers $\{u_{i}\}_{i=1}^{|V|}$ and $\{v_{j}\}_{i=1}^{|V|}$ such that $|\mathcal{M}_{ij}-u_{i}v_{j}|\leq 3/2-0.001|V|^{-6}$ for all $i,j$.

Proof.

Assume $(G,D)$ admits a valid coloring $(W,B)$ as in Problem 3, and let $G^{\prime}$ be the directed graph obtained after the transformation in item (ii) of Problem 3. Consider the matrix $N$ defined by $N_{ij}=\mathcal{M}_{ij}$ if $\{i,j\}$ is monochromatic and $N_{ij}=-\mathcal{M}_{ij}$ otherwise. The matrix $N$ has $2$’s on the diagonal, $-1$’s at every position $(i,j)$ which is an edge of $G^{\prime}$, and zeros and ones everywhere else. Since $G^{\prime}$ has no cycle, by Observation 1, there is a permutation matrix $C$ (corresponding to an ordering) such that all the $-1$’s are located below the main diagonal of the matrix $C^{-1}NC$. Since the permutations of rows and columns and multiplications of them by $-1$ do not change our ability or inability to approximate a matrix, it is sufficient to prove the existence of real numbers $\{u_{i}\}_{i=1}^{|V|}$ and $\{v_{j}\}_{i=1}^{|V|}$ such that $|M_{ij}-u_{i}v_{j}|\leq 3/2-0.001|V|^{-6}$ for any $n\times n$ matrix $M$ with the $2$’s on the diagonal, $-1$’s at some positions below the diagonal, and zeros and ones everywhere else. Towards this end, one can check that it suffices to set

[TABLE]

where $\varepsilon=0.1|V|^{-4}$.

Conversely, assume that the numbers $\{u_{i}\}_{i=1}^{|V|}$ and $\{v_{j}\}_{i=1}^{|V|}$ are such that $|\mathcal{M}_{ij}-u_{i}v_{j}|\leq 3/2-0.001|V|^{-6}<3/2$ for all $i,j$. Since $\mathcal{M}_{ii}=2$ for all $i$, we have $|u_{i}v_{i}-2|<3/2$ which implies that $u_{i}$ and $v_{i}$ are non-zero and have the same sign. A relabeling of indices does not change the properties we discuss, so we can assume $|v_{1}|\leq\ldots\leq|v_{n}|$. We define $W$ (resp. $B$) as the set of all $i$ such that $u_{i}>0$ (resp. $u_{i}<0$), and we are going to prove that $(W,B)$ is a valid coloring of $G$ as in Problem 3. We define the matrix $N$ by multiplying the rows and columns of $\mathcal{M}$ with indices in $B$ by $-1$, and define $u^{\prime}$ and $v^{\prime}$ by multiplying the entries in $u$ and $v$, respectively, with indices in $B$ by $-1$. We have $\left|N_{ij}-u^{\prime}_{i}v^{\prime}_{j}\right|<3/2$ and $u^{\prime}\geq 0$ and $v^{\prime}\geq 0$. For $j>i$, $v^{\prime}_{j}\geq v^{\prime}_{i}$ which implies $u^{\prime}_{i}v^{\prime}_{j}\geq u^{\prime}_{i}v^{\prime}_{i}>1/2$, so $N_{ij}\neq-1$ if $j>i$. In other words, all the $-1$’s of $N$ are located below the main diagonal, and we have $N_{ab}=N_{ba}=1$ for all the positions $(a,b)$ in $D$. This means that the item (i) of Problem 3 is satisfied, and the graph as in item (ii) has indeed no cycle because the edges $(a,b)$ of this graph correspond to the positions of the $-1$’s in $N$ which are all located below the main diagonal. ∎

Now we can determine the complexity status of D-$ell_{\infty}$-R1A by proving that it is NP-complete.

Theorem 3**.**

The D-$ell_{\infty}$-R1A problem is NP-complete.

Proof.

The function $(G,D)\to(\mathcal{M},3/2-0.001|V|^{-6})$ can be computed in polynomial time, and Theorem 2 proves that it is a reduction from Problem 3 to D-$ell_{\infty}$-R1A. Therefore, D-$ell_{\infty}$-R1A is NP-hard by Lemma 3. The membership of D-$ell_{\infty}$-R1A in NP is stated in Theorem 1. ∎

Example 2**.**

For the matrix

[TABLE]

there exist a permutation and a sign flip of the rows and columns so that the negative entries are below the main diagonal: $N_{ij}=M_{\pi_{i}\pi_{j}}s_{\pi_{i}}s_{\pi_{j}}$ with $s=(1,-1,1,1,1)$ and $\pi=(5,2,1,4,3)$, with

[TABLE]

We have $\min_{u,v}||M-uv^{T}||_{\infty}=1.3456<3/2$. For the matrix

[TABLE]

there is no such sign flip and permutation, and $\min_{u,v}||M-uv^{T}||_{\infty}=3/2$.

4 Heuristic algorithm and application to the recovery of quantized low-rank matrices

This goal of this section is to describe a simple heuristic algorithm for $\ell_{\infty}$ LRA and apply it for the recovery of quantized low-rank matrices. It will allow us to get some more insight on this problem. The algorithm is available from https://sites.google.com/site/nicolasgillis/ and allows the interested reader to tackle $\ell_{\infty}$ LRA (in particular the examples presented in this paper can be run directly). In this section, all tests are performed using Matlab on a laptop Intel dual CORE i5-3210M CPU @2.5GHz 6Go RAM.

4.1 Block coordinate descent method

A popular approach in optimization is block coordinate descent (BCD): fix a subset of the variables and optimize over the other variables; see [22] for a recent survey. An crucial aspect of BCD is to make the subproblem easy (and fast) to solve. For $\ell_{\infty}$ LRA, a judicious choice is to optimize alternatively over the columns of $U$ and the rows of $V$; see Algorithm 1. In fact, the corresponding subproblems are convex and separable, that is, each entry in a column of $U$ (resp. in a row of $V$) can be optimized independently of the other entries in the same column (resp. row); this is described in the next section.

To initialize Algorithm 1, we use the optimal solution of $\ell_{2}$ LRA. It would be an interesting direction of research to use more sophisticated initialization strategies such as the approximation algorithm proposed in [2] that can be refined by Algorithm 1.

4.2 Secant method for the subproblem

Let us focus on the rank-one subproblem in $v$ (by symmetry, the subproblem in $u$ can be solved in the same way). It can be decoupled into $n$ problems in one variable: for $1\leq j\leq n$, we need to solve

[TABLE]

The optimal solution is not necessarily unique. Non-uniqueness may happen when $u_{i}=0$ and $|M_{ij}|$ is large for some $i$, while uniqueness is guaranteed if $u_{i}\neq 0$ for all $i$ because the objective function is piece-wise linear with nonzero slopes. To make the problem well posed, it makes sense to consider

[TABLE]

with a unique solution which is also optimal for (4.1). (Note that if $u=0$, any $v_{j}$ is optimal.) Let us focus w.l.o.g. on the case $u\geq 0$ by flipping the signs of the rows of $M$ accordingly. The global minima is the intersection of two linear functions of the form $f^{i}_{\pm}(v_{j})=\pm(M_{ij}-u_{i}v_{j})$ ($1\leq i\leq m$), one with negative slope and one with positive slope. Therefore, the optimal solution $v_{j}^{*}$ of (4.2) satisfies

[TABLE]

Hence solving (4.2) can be done by identifying the two indices $i_{1}$ and $i_{2}$ corresponding to the optimal solution. This could be performed by inspection since there are $\frac{m(m-1)}{2}$ possible pairs. A more efficient approach is described in the following. Let us define

[TABLE]

We have $v_{l}=\frac{M_{i_{l}j}}{u_{i_{l}}}\leq v_{j}^{*}\leq\frac{M_{i_{u}j}}{u_{i_{u}}}=v_{u}$. Since the objective function is convex, it is rather straightforward to implement the following secant method:

Initialize $(i_{1},i_{2})=(i_{l},i_{u})$. 2. 2.

Intersect the two lines corresponding to the indices $i_{1}$ (with negative slope) and $i_{2}$ (with positive slope) to obtain the point $v_{c}=\frac{M_{i_{1}j}+M_{i_{2}j}}{u_{i_{1}}+u_{i_{2}}}$. 3. 3.

Compute the objective function of (4.2) in $v_{j}=v_{c}$ and identify the index $i_{a}$ that is active, that is, the index $i_{a}$ such that $u_{i_{a}}\neq 0$ and $|M_{i_{a}j}-u_{i_{a}}v_{j}|=\max_{\{i|u_{i}\neq 0\}}|M_{ij}-u_{i}v_{j}|$. If the slope in $i_{a}$ is positive, replace $i_{2}$ by $i_{a}$; otherwise replace $i_{1}$ by $i_{a}$. If the two indices $(i_{1},i_{2})$ are active together, the algorithm has converged: return $v_{j}^{*}=\frac{M_{i_{1}j}+M_{i_{2}j}}{u_{i_{1}}+u_{i_{2}}}$; otherwise, return to 2.

It turns out that this secant method performs surprisingly well in the sense that it needs a very small number of iterations to terminate. Let us illustrate this on randomly generated instances.

Example 3** (Numerical experiment on the secant method for (4.1)).**

We have run the above secant method to solve $10^{4}$ problems of the form (4.1) for different values of $m$, generating each entry of $M$ and $u$ using the normal distribution $N(0,1)$. Table 1 reports the distribution of the number of iterations needed to solve (4.1) for the $10^{4}$ problems, along with the total computational time to solve them.

*We observe that the secant method requires in average 5 iterations to terminate while 9 are necessary in the worst case. *

Computational cost

Each iteration of the secant method requires $O(m)$ operations hence solving (4.2) requires $O(mK)$ operations where $K$ is the number of iterations performed by the secant method. In practice, as illustrates by Table 1, $K$ will be rather small (in all cases we have generated randomly, $K$ is smaller than 9). In the worst case, the secant could potentially have to go through all indices with $K=O(m)$ and require $O(m^{2}n)$ operations. Note that sorting the values $\frac{M_{ij}}{u_{i}}$ ($1\leq i\leq m$) in $O(m\log(m))$ operations would allow to perform a bisection in $O(m\log(m))$ operations (getting rid of about $m/2$ indices per iteration). Finally, the total computational cost of Algorithm 1 is $O(mnrK)$ per iteration.

Remark 2** ($\ell_{\infty}$ LRA with nonnegativity constraints).**

Algorithm 1 can easily be adapted to incorporate nonnegativity constraints on $U$ and $V$. In fact, the optimal solution of (4.1) with the constraint $v_{j}\geq 0$ is $\max(v_{j}^{*},0)$ where $v_{j}^{*}$ is the optimal solution of the unconstrained problem (4.1).

4.3 Recovery of quantized low-rank matrices

We apply in this section Algorithm 1 to recover a quantized low-rank matrix. First, let us start with a toy example.

Example 4**.**

Let us generate a simple example with $m=8$, $n=5$ and $r=3$ where $M=UV$ with each entry of $U$ and $V$ generated using the normal distribution $N(0,1)$. We obtain (with two digits of accuracy)

[TABLE]

and a quantization (here we simply use its rounding to the nearest integer)

[TABLE]

which has rank 5, and with $||M-M_{q}||_{\infty}=0.498$ (entry on second row, third column). The optimal solution $X_{\text{svd}}$ of rank-3 $\ell_{2}$ LRA gives $||M_{q}-X_{\text{svd}}||_{\infty}=0.57$ while Algorithm 1 provides a solution $X^{*}$ with $||M_{q}-X^{*}||_{\infty}=0.39$. Note that, for this problem, the set

[TABLE]

is rather large and does not only contain a small neighborhood around the matrix $M$. It would be an interesting direction for further research to identify conditions for this problem to be well posed (e.g., sufficiently many entries of $M-M_{q}$ close to $\pm 1/2$, or using additional constraints on $X$).

Let us construct instances exactly as in Example 4 except that we use $m=n=200$ and $r=1,2,5,10,20$. Note that the advantage of using quantized low-rank matrices is that we know there there exists a solution with error smaller than 1/2. We stop the execution of Algorithm 1 only when 1000 iterations are performed or when the relative error between two iterates is smaller than $10^{-6}$, that is, when $e_{t}-e_{t+1}\leq 10^{-6}||M_{q}||_{\infty}$, where $e_{t}$ is the error at iteration $t$. Table 2 provides the smallest, average and largest value of the error of Algorithm 1 (second column), the number of solutions found with error smaller than 0.5 (third column), the smallest, average and largest number of iterations needed to converge (fourth column), the average computational time (fifth column), and the smallest, average and largest value of the $\ell_{\infty}$ error for $\ell_{2}$ LRA which we used as an initialization for Algorithm 1 (last column).

We observe the following:

• •

In all cases, Algorithm 1 is able to significantly improve the initial solution computed with $\ell_{2}$ LRA (second column vs. last column of Table 2).

• •

In many cases, Algorithm 1 converges in a relative few number of iterations (sometimes in 3 iterations); see the fourth column of Table 2. We believe the reason is that the objective function landscape is rather peaky hence it is more likely for the algorithm to terminate rapidly.

• •

For $r=1$, Algorithm 1 is always able to recover a solution with error smaller than 1/2. This is not surprising since the problem is not difficult: in fact, the graph $G_{b}(M_{q},0.5)$ contains a single connected component since most entries of $|M_{q}|$ are larger than 0.5 (see Theorem 1). We observed in practice that, in this case, Algorithm 1 is always able to converge to an optimal solution. Optimality can be verified by checking that the answer to D-$\ell_{\infty}$-R1A($M_{q}$,$f^{*}-\epsilon$) is NO, where $f^{*}$ is the objective function value of Algorithm 1 at convergence and $\epsilon$ is a small positive constant. (We have included this function in the code available online.)

• •

When $r>1$, the problem becomes more difficult and Algorithm 1 is not able to identify a solution with error smaller than 1/2 (it is never able to do it for $r=2,5,10$). This leads to an interesting question: is $\ell_{\infty}$ LRA hard for $r>1$ (say $r=2$) even when $G_{b}(M_{q},0.5)$ contains a unique connected component (or when $M\geq 0$)? Table 2 suggests that it is the case. At least, we observed that there are many local minima: for example, for $r=2$, using 100 random initializations does not allow in general to obtain a solution with error smaller than 1/2, while the solution $M$ cannot be improved; see Table 3. Note also that the solutions obtained with random initializations have error significantly larger than with $\ell_{2}$-LRA initialization.

• •

Surprisingly, when $r=20$, Algorithm 1 is able in most cases to recover a solution with error smaller than 1/2. We believe the reason is that the number of degrees of freedom is large hence the optimal solution has error smaller than 0.5. This is confirmed by the results in Table 3 where we have initialized Algorithm 1 using the original rank-$r$ matrix $M$ (hence the initial error is close to 0.5). We see that, quite naturally, it is able to identify better solutions than when initialized with the solution of $\ell_{2}$ LRA.

5 Conclusion

In this paper, we have analyzed the component-wise $\ell_{\infty}$ low-rank matrix approximation problem. We proved that, even in the rank-one case, the decision version of this problem is NP-complete using a reduction from ‘NOT ALL EQUAL 3SAT’ (Theorem 3). However, in the rank-one case when $M\geq 0$ or when $G_{b}(M,k)$ contains a few number of connected components, the problem can be solved in polynomial time (Theorem 1). We then described a simple block coordinate descent method and applied it for the recovery of quantized low-rank matrices. We observed that, as expected, the algorithm is able to recover an optimal solution when $r=1$. However, as soon as $r>1$, the problem becomes more difficult even when $G_{b}(M,k)$ contains a single connected component.

Acknowledgment

We are grateful to Laurent Jacques for pointing us out to the quantized low-rank matrix problem and for insightful discussions. Part of this paper was written at the University of Mons during a visit of the second author who is grateful for the invitation and for being introduced to this problem.