Approximating $q\to p$ Norms of Non-Negative Matrices in Nearly-Linear Time

We give the first input-sparsity time algorithm to approximate the ${\mathrm{\ell }}_{q\to p}$-norm of a non-negative matrix when $q\ge p\ge 1$. The result is stated in the following theorem.

While the statement of Theorem 1 concerns approximation of a decision problem, standard techniques can transform an approximation decision algorithm into an approximation algorithm. We refer the reader to Appendix A for the details which we provide for completeness.

In the context of tree-based oblivious routing schemes, the competitive ratio can be written directly as an ${\mathrm{\ell }}_{p\to p}$-norm expression. If $𝑴$ is the $\left|E\right|\times \left|E\right|$ matrix that specifies how unit demands on edges are routed, then the competitive ratio is ${\max }_{{\Vert 𝒙\Vert }_p\le 1}{\Vert 𝑴𝒙\Vert }_p$. [4] proved it is sufficient to approximate ${\Vert 𝑴\Vert }_{p\to p}$ a polynomial number of times to compute a $O(\log n)$ competitive tree-based oblivious routing. Hence, we can use Theorem 1 for the case $p=q$, and $𝑴\in {\mathrm{ℝ}}^{m\times m}$ to improve the running time of [4] by a factor of $O(m^3)$. Moreover, if both the load function and the aggregation function are known monotone norms, we provide the first algorithm to compute an optimal oblivious linear routing. Notice that this problem is slightly different from the one solved by [4] as we do need knowledge of the load function. Hence, the following theorem.

We seek to maximize $g:𝒙\mapsto {\Vert 𝑨𝒙\Vert }_p^q/{\Vert 𝒙\Vert }_q^q$. Our approach is partly inspired from the long line of work on solvers for positive linear programs [17, 24, 1, 18], as well as certain non-standard instantiations of the same framework in the context of regression problems [8]. Those algorithms are width-independent, meaning their running time is at most linear in $\text{nnz}(𝑨)+n+m$ (at a cost of a larger dependency on the precision $\epsilon$). The width of an instance measures the instance’s diameter in the appropriate norm. For example, for linear programs, the width is typically defined as the maximum absolute value that any of its linear constraint functions can attain over the search space. By design, those algorithms are different from classical multiplicative weight update algorithms as their running time does not depend on the width of the problem.

At each iteration, we scale by $(1+\alpha )$ some coordinates of our iterate to ensure large progress is made. We say that a coordinate is hit at time $t$ if we scale it between iterations $t$ and $t+1$. We do not directly measure the progress on $g$, instead, we separate the progress obtained in the denominator from the progress obtained in the numerator. At each iteration, assuming $𝑨$ is non-negative, we want to ensure the ratio between those progresses is larger than ${(1-O(\epsilon ))}^q{\Vert 𝑨\Vert }_{q\to p}^q$ (7 with $V={\Vert 𝑨\Vert }_{q\to p}$). This ensures after $T$ iterations that $g({𝒙}^{(T)})$ is a $(1-O(1/{\Vert {𝒙}^{(T)}\Vert }_q^q))$ approximation of ${\Vert 𝑨\Vert }_{q\to p}^q$. To do so, we consider a potential function $𝚽:{\mathrm{ℝ}}^n\to {\mathrm{ℝ}}^n$ (8) and hit the coordinates whose potential is above a threshold. $𝑨$ being non-negative is crucial as it keeps the property that the norm $𝒙\mapsto {\Vert 𝑨𝒙\Vert }_p$ is monotone.

The algorithm converges once ${\Vert {𝒙}^{(T)}\Vert }_q$ is large. If we were to hit only one coordinate of ${𝒙}^{(t)}$ at iteration $t$ (say the coordinate with the largest potential), the running time would depend on the ambient dimension $n$ which would not be a width-independent algorithm. To overcome this issue, we transform our problem into an approximate decision one. Given $V\ge 0$, a guess on ${\Vert 𝑨\Vert }_{q\to p}$, we either want to find $𝒙$ such that $g(𝒙)\ge {(1-\epsilon )}^q{V}^q$ or certify that $V>{\Vert 𝑨\Vert }_{q\to p}$. We then change the threshold for $𝚽$ accordingly. The certification comes from the fact that at each iteration, at least one coordinate of $𝚽(𝒙)$ should be above ${\Vert 𝑨\Vert }_{q\to p}^q$ (10). Assuming $V\le {\Vert 𝑨\Vert }_{q\to p}$, there should always be a coordinate with potential larger than $V^q$. This modification still does not give a width-independent algorithm, however, $𝚽$ behaves “nicely” under multiplicative scaling which allows us to upper bound coordinate-wise the value of $𝚽({𝒙}^{(t+1)})/𝚽({𝒙}^{(t)})$ (Lemma 11).

We use this upper bound to prevent new coordinates from having potential larger than $V^q$. This allows us to guarantee that at iteration $t$, a coordinate with potential larger than $V^q$ also had potential larger than $V^q$ at iterations $0,1,\mathrm{\dots },t-1$. Since its potential was always larger than $V^q$, it was always hit and so it is large. In order to prevent coordinates from having potential larger than $V^q$, we need to also hit coordinates with potential slightly smaller than $V^q$ which is not an issue since we aim for an approximation. This gives a first width-independent algorithm to approximate the ${\mathrm{\ell }}_{q\to p}$-norm of a non-negative matrix in time $\tilde{O}(\text{nnz}(𝑨)/q\epsilon )$ or in parallel time $\tilde{O}(1/q\epsilon )$ (Corollary 14).

The algorithm can also use preconditioning techniques for ${\mathrm{\ell }}_p$ subspace embedding. We first compute a matrix ${𝑨}^{\prime }$ such that ${\Vert 𝑨𝒙\Vert }_p=(1\pm \epsilon ){\Vert {𝑨}^{\prime }𝒙\Vert }_p$ and ${𝑨}^{\prime }$ has fewer nonzero entries than $𝑨$. Lewis weights sampling introduced by Cohen and Peng [7] and improved for $p>2$ in [23] allows to reduce the number of rows of $𝑨\in {\mathrm{ℝ}}^{m\times n}$ to $\tilde{O}(n/{\epsilon }^2)$ when $p\le 2$. This preconditioning technique, paired with the fact that ${\Vert 𝑨\Vert }_{q\to p}={\Vert {𝑨}^T\Vert }_{p^{\ast }\to q^{\ast }}$ where $1/q+1/q^{\ast }=1/p+1/p^{\ast }=1$, allows better running time when $m\gg n$ or $n\gg m$. This result is showed in Corollary 15.

The algorithm is linked with the Perron-Frobenius theorem, in the sense that for $q=p=2$, if $𝑨$ is positive then the maximum of $𝒙\mapsto {\Vert 𝑨𝒙\Vert }_2/{\Vert 𝒙\Vert }_2$ is the spectral norm of $𝑨$. For positive matrices, 10 is a generalized version of the Collatz-Wielandt formula which states that ${\min }_{𝒙>\mathrm{𝟎}}{\max }_i\frac{⟨{𝑨}_{\bullet ,i},𝑨𝒙⟩}{{𝒙}_i}={\Vert 𝑨\Vert }_{2\to 2}^2$. We generalize this formula to ${\mathrm{\ell }}_p$ norms, and obtain ${\min }_{𝒙>\mathrm{𝟎}}{\max }_i\frac{⟨{𝑨}_{\bullet ,i},{(𝑨𝒙)}^{p-1}⟩}{{𝒙}_i^{p-1}}={\Vert 𝑨\Vert }_{p\to p}^p$. When $p\ne q$, we just add a scaling factor that depends on ${\Vert 𝒙\Vert }_q$ to find the equality.

We then show that an application of our algorithm is to compute a $O(\log n)$-competitive linear oblivious routing with a better running time than [4]. When the load function is a known monotone norm, we show that we can reduce the problem to finding a saddle point of a bi-linear function.

Computing the ${\mathrm{\ell }}_{q\to p}$ norm on general matrices is hard. For arbitrary matrices, there are three known easy cases: ${\Vert 𝑨\Vert }_{2\to 2}$, ${\Vert 𝑨\Vert }_{q\to \mathrm{\infty }}$ and ${\Vert 𝑨\Vert }_{1\to p}$ (see [22]).

On the hardness front, for arbitrary matrices, computing ${\Vert 𝑨\Vert }_{q\to p}$ when is NP-hard [22]. For $q=p\ne 1,2,\mathrm{\infty }$, even if $p$ and the matrix are assumed to have rational entries, it is NP-hard to compute ${\Vert 𝑨\Vert }_{p\to p}$ [12].

If we aim to find an approximation, for arbitrary matrices the problem is also hard to solve to arbitrary precision. For $q\ge p>2$ (and $2>q\ge p>1$), assuming $\mathrm{NP}\notin \mathrm{DTIME}(n^{\mathrm{poly}\log (n)})$, the problem cannot be approximated to a factor $2^{{(\log n)}^{1-\epsilon }}$, for any constant $\epsilon >0$ [4]. [3] showed that under the Small-Set Expansion hypothesis, it is NP-hard to approximate ${\Vert 𝑨\Vert }_{2\to p}$ for arbitrary matrices $𝑨$ when $p\ge 4$. Finally, when (and ), if $\mathrm{NP}\notin \mathrm{BPTIME}(2^{{(\log n)}^{O(1)}})$, it is NP-hard to approximate ${\Vert 𝑨\Vert }_{q\to p}$ to a factor $2^{{(\log n)}^{1-\epsilon }}$ for any constant $\epsilon >0$ [5].

For non-negative matrices, [6] used a power iteration-type algorithm to approximate $p$-norm of non-negative matrices without any bounds on the time of convergence. When $q\ge p\ge 1$, [4] provided the analysis and proved that the power iteration-type algorithm from [6] converges in time $\tilde{O}(\frac{n{(n+m)}^2}{\epsilon }\text{nnz}(𝑨))$. Under the same settings, [22] demonstrated that the problem is equivalent to maximizing a concave function ($𝒙\mapsto {\Vert 𝑨{𝒙}^{\frac{1}{q}}\Vert }_p$) over the unit simplex. When , no hardness result is known for non-negative matrices.

The non-negative assumption may be necessary to make the problem easy. Indeed, it is possible to relax MAXCUT into approximating the ${\mathrm{\ell }}_{\mathrm{\infty }\to 1}$ norm of a symmetric diagonally dominant matrix (see [22]). Since MAXCUT is known to be NP-Had to approximate better than $16/17\approx 0.941$ [13], and UG-Hard to approximate better than $\approx 1/0.878$ [16], this implies that there is no PTAS for the ${\mathrm{\ell }}_{\mathrm{\infty }\to 1}$ problem, even for symmetric diagonally dominant matrices.

Hence, approximating ${\Vert 𝑨\Vert }_{\mathrm{\infty }\to 1}$ is NP-hard even when $𝑨$ is positive semi-definite and the precision sought is fixed. However, we can easily extend the results on non-negative matrices to matrices $𝑨$ such that there exist two diagonal matrices $𝑳\in {\mathrm{ℝ}}^{m\times m}$ and $𝑹\in {\mathrm{ℝ}}^{n\times n}$ with $\pm 1$ entries on their diagonals such that $\mathrm{𝑳𝑨𝑹}$ is non-negative [22]. Indeed, given $𝑳,𝑹$, we have that for any $𝒙$, ${\Vert 𝑨𝒙\Vert }_p={\Vert \mathrm{𝑳𝑨}𝒙\Vert }_p$. Hence, we can apply the algorithm on $\mathrm{𝑳𝑨𝑹}$ and consider $𝑹𝒙$ where $𝒙$ is a $(1-\epsilon )$-approximation of ${\Vert \mathrm{𝑳𝑨𝑹}\Vert }_{q\to p}$.

An immediate application of approximating the ${\mathrm{\ell }}_p$ norm of non-negative matrices is computing oblivious routings, per [9]. The main difficulty posed by their approach was efficiently certifying whether a non-negative matrix implicitly defined by the underlying graph has large or small induced ${\mathrm{\ell }}_p$ norm, for $p>2$. The polynomial time algorithm of [4] could be used together with the [9] framework, at the expense of paying a large overall running time. Notably, for the case of $p=2$, approximating the induced norm can be solved efficiently via the power method; this special case has been leveraged by [10] to achieve ${\mathrm{\ell }}_{\mathrm{\infty }}$ oblivious routings running in sub-quadratic time.

We use lowercase bold letters for vectors and uppercase for matrices. For a vector $𝒛$ (resp. a matrix $𝑨$), ${𝒛}^T$ (resp. ${𝑨}^T$) denotes the transpose of $𝒛$ (resp. of $𝑨$). For an integer $i$, ${𝑨}_i$ denotes the $i$^th row of $𝑨$, and ${𝑨}_{\bullet ,i}$ the $i$^th column. We say that a vector or a matrix is non-negative if all its entries are non-negative. $⟨.,.⟩$ is the usual inner product and ${⟨.,.⟩}_F$ denotes the Frobenius inner product. For $p\ge 1$ and a vector $𝒛$, we denote by ${\Vert 𝒛\Vert }_p$ the ${\mathrm{\ell }}_p$-norm of $𝒛$ defined as

For $q,p\ge 1$ and a matrix $𝑨$, we denote by ${\Vert 𝑨\Vert }_{q\to p}$ the ${\mathrm{\ell }}_{q\to p}$-norm of $𝑨$ defined as

$\mathrm{𝟎}$ (resp. $\mathrm{𝟏}$) denotes the vector with all entries equal to $0$ (resp. $1$). For a vector $𝒛$ and a real $q$, ${𝒛}^q$ is the vector where the exponent is applied coordinate wise. For a matrix $𝑨$, $\text{nnz}(𝑨)$ denotes the number of nonzero elements in $𝑨$. We denote by ${\mathrm{\Delta }}_n$ the non-negative simplex ${\mathrm{\Delta }}_n:=\{𝒛\in {\mathrm{ℝ}}_{\ge 0}^n:{\sum }_i{𝒛}_i=1\}$.

We would like to emphasize that the problem of maximizing $f(𝒙)={\Vert 𝑨𝒙\Vert }_p$ on the unit ${\mathrm{\ell }}_q$ ball corresponds to maximizing a convex function over a convex domain which can not be solved using classical convex optimization tools. Nevertheless, when $𝑨$ is non-negative, $q\ge p\ge 1$, we can show that by composing $f$ with a simple concave mapping from the unit simplex ${\mathrm{\Delta }}_n$ to the unit ${\mathrm{\ell }}_q$ sphere we obtain a concave function that we would like to maximize over a convex domain [22]. To intuit, one should remark that the maximum of $f$ is achieved by a vector in the positive orthant of the unit ${\mathrm{\ell }}_q$ sphere which we call $S_q$. Moreover, $S_q$ is not a convex set. In fact, for two different vectors $𝒙,𝒚\in S_q$, and , we have , hence we are “losing” norm by taking the convex combination of two points. That means that even if $f$ is convex, we have the potential to “increase” the value of $f$ between any two points of $S_q$ by taking a path in $S_q$. This new function has better chance to be concave than $f$ since its value between two points of $S_q$ is larger. One can remark that if $𝒖\in {\mathrm{\Delta }}_n$, then ${𝒖}^{\frac{1}{q}}\in S_q$. The mapping $g:𝒖\mapsto {𝒖}^{\frac{1}{q}}$ describes such a path since for any $t\in [0,1]$ and $𝒖,𝒗\in {\mathrm{\Delta }}_n$, we have that $tg(𝒖)+(1-t)g(𝒗)\in S_q$. As expected, $g$ is entry-wise concave and calculus gives that as long as $q\ge p$, $𝒖\mapsto f(g(𝒖))={\Vert 𝑨{𝒖}^{\frac{1}{q}}\Vert }_p$ is concave. We give the proof in Subsection D.1.

This shows that a simple re-parametrization of the problem turns it into a concave maximization problem, which can be efficiently solved using a cutting plane method [14], leading to a running time of $\tilde{O}(n\cdot \text{nnz}(𝑨)+n^3)$. However, as cubic runtime may turn out to be prohibitive, we will focus on more efficient algorithms at the expense of a linear dependence in the error tolerance of the provided solution.

The previous state-of-the-art algorithm was provided in [4]. Let

In order to compute ${\Vert 𝑨\Vert }_{q\to p}$, one needs to maximize $f$ over ${\mathrm{ℝ}}^n$. If we assume $𝑨$ is non-negative, then [4] provide an algorithm with running time $\tilde{O}(\frac{n{(n+m)}^2}{\epsilon }\text{nnz}(𝑨))$. To find the maximum of $f$, we can compute its gradient $\nabla f$. We have

At the optimum, $\nabla f=\mathrm{𝟎}$ which can be written as

Since $𝑨$ is non-negative, we know there is a maximum $𝒙$ with non-negative coordinates. [4] then define the operator $S:{\mathrm{ℝ}}_{\ge 0}^n\to {\mathrm{ℝ}}_{\ge 0}^n$ such that $S(𝒙)={\left({𝑨}^T{(𝑨𝒙)}^{p-1}\right)}^{\frac{1}{q-1}}$. It is easy to show that if $𝒙$ is such that $S(𝒙)\propto 𝒙$, then $𝒙$ is a critical point and ${\Vert 𝑨\Vert }_{q\to p}=\frac{{\Vert 𝑨𝒙\Vert }_p}{{\Vert 𝒙\Vert }_q}$. Consider the two potentials $m(𝒙):={\min }_{i}S{(𝒙)}_i/{𝒙}_i$ and $M(𝒙)={\max }_{i}S{(𝒙)}_i/{𝒙}_i$, [4] showed the following lemma.

At optimum, $m{({𝒙}^{\ast })}^{q-1}=M{({𝒙}^{\ast })}^{q-1}={\Vert 𝑨{𝒙}^{\ast }\Vert }_p^p/{\Vert {𝒙}^{\ast }\Vert }_q^q$. Hence, the goal is to reduce the ratio $M(𝒙)/m(𝒙)$. [4] showed that iteratively applying $S$ such that ${𝒙}^{(t+1)}=S({𝒙}^{(t)})$ outputs ${𝒙}^{(T)}$ a $(1-\epsilon )$-approximation when $T=\tilde{O}(\frac{n{(n+m)}^2}{\epsilon })$. Hence, the total running time is $\tilde{O}(\frac{n{(n+m)}^2}{\epsilon }\text{nnz}(𝑨))$. In this paper, we provide a much simpler algorithm with running time $\tilde{O}(\text{nnz}(𝑨)/q\epsilon )$.

Cohen and Peng introduced in [7] a fast algorithm to compute a matrix ${𝑨}^{\prime }$ such that ${\Vert 𝑨𝒙\Vert }_p{\approx }_{1+\epsilon }{\Vert {𝑨}^{\prime }𝒙\Vert }_p$ for any vector $𝒙$. The algorithm computes ${𝑨}^{\prime }$, a matrix containing few rescaled rows of $𝑨$. The algorithm is of particular interest in our situation when $𝑨\in {\mathrm{ℝ}}^{m\times n}$ is such that $m\gg n$. For $1\le p\le 2$, $\tilde{O}(\frac{n}{{\epsilon }^2})$ rows are sufficient. Moreover, the cost of such an algorithm has only time complexity $\tilde{O}(\text{nnz}(𝑨)+n^{\omega })$ where $\omega$ is the time complexity of matrix multiplication. Formally, this gives the following theorem.

Since the algorithm just samples and rescales rows of $𝑨$, if $𝑨$ is non-negative, then so is ${𝑨}^{\prime }$. Hence, it makes Lewis weights sampling a possible preconditioning of the input.

We want to solve the following approximate decision problem, which we show is sufficient to solve the approximate optimization problem in Appendix A.

We seek to maximize a function $f$ of the form $f:𝒙\mapsto \frac{f_n(𝒙)}{f_d(𝒙)}$ where both $f_n$ and $f_d$ are convex functions. Hence, we have

If we take $f_d:𝒙\mapsto {\Vert 𝒙\Vert }_q$, we have $\nabla f_d(𝒙+𝜹)={\left(\frac{𝒙+𝜹}{{\Vert 𝒙+𝜹\Vert }_q}\right)}^{q-1}$ which is not optimal to upper bound by $\nabla f(𝒙)$. Instead, we consider $g_n:𝒙\mapsto {\Vert 𝑨𝒙\Vert }_p^q$, $g_d:𝒙\mapsto {\Vert 𝒙\Vert }_q^q$ and we want to maximize $g:𝒙\mapsto {\left(\frac{{\Vert 𝑨𝒙\Vert }_p}{{\Vert 𝒙\Vert }_q}\right)}^q$. We now obtain a tighter bound

This bound will be used to design an iterate that satisfies the following lemma. We prove 7 in Subsection E.1.

With scaling updates, i.e. ${𝒙}^{(t+1)}={𝒙}^{(t)}+{𝜹}^{(t)}$ where ${𝜹}_i^{(t)}\in \{0,\alpha {𝒙}_i^{(t)}\}$ we obtain the following inequality.

To simplify the notation, we define a vector of potentials $𝚽(𝒙)\in {\mathrm{ℝ}}^n$. Those potentials have nice properties, for instance, the direction of the gradient of $g$ can be extracted from the values of the potentials. Indeed, they are made such that $\nabla g(𝒙)=\frac{\nabla g_d(𝒙)}{g_d(𝒙)}\circ \left(𝚽(𝒙)-g(𝒙)\mathrm{𝟏}\right)$. That means when the potential is larger than $g(𝒙)$, the gradient is positive.

The algorithm does not follow the traditional gradient descent framework, however each iteration still follows the direction of the gradient. We want our iterate to satisfy the following corollary of 7.

For the iterate to satisfy the invariant in 9, there must always be a potential larger than $V^q$. This is guaranteed by the following lemma from [4].

As long as $V\le {\Vert 𝑨\Vert }_{q\to p}$, there always is a coordinate with potential above $V^q$. Setting $\alpha =\epsilon /2$, at time $t$ we scale coordinates with potential larger than $V^q$ so that our iterate satisfies Invariant (2). At each iteration, a different coordinate may be scaled which means ${\Vert {𝒙}^{(T)}\Vert }_q$ increases quite slowly. We only have the following lower bound ${\Vert {𝒙}^{(T)}\Vert }_q^q\ge {(1+\alpha )}^{T/n}$, which requires $T=\tilde{O}(n/q\epsilon )$ iterations are required to solve the approximate decision problem.

In order to achieve an efficient algorithm, it is necessary to improve the number of iterations by a factor of $n$. The challenge we had to deal before was that the ${\mathrm{\ell }}_q$ norm of ${𝒙}^{(t)}$ was not increasing rapidly enough. To overcome this issue, we will modify the algorithm to ensure there is a coordinate $k$ such that ${𝒙}_k^{(t)}$ is large. Notice that 10 guarantees that there is at least one coordinate with potential larger than $V^q$. The new invariant we would like to implement is that the set of coordinates with potential larger than $V^q$ does not take any new elements. This requires to also scale coordinates with a potential slightly smaller than $V^q$ which does not prevent satisfying the first Invariant (2). In order to do so, we need to analyze how the potentials evolve between two iterations.