Efficient nuclear norm approximation via the randomized UTV algorithm

TL;DR

This paper presents a modified randUTV algorithm optimized for efficiently approximating the singular values of a matrix, enabling accurate nuclear norm computation with minimal additional cost.

Contribution

It introduces a variant of randUTV tailored for singular value approximation and nuclear norm estimation, improving efficiency over the original full factorization approach.

Findings

Effective approximation of singular values using the modified randUTV

Accurate nuclear norm estimation with negligible additional cost

Enhanced computational efficiency for Schatten-$p$ norms

Abstract

The recently introduced algorithm randUTV provides a highly efficient technique for computing accurate approximations to all the singular values of a given matrix $A$. The original version of randUTV was designed to compute a full factorization of the matrix in the form $A = UTV^*$ where $U$ and $V$ are orthogonal matrices, and $T$ is upper triangular. The estimates to the singular values of $A$ appear along the diagonal of $T$. This manuscript describes how the randUTV algorithm can be modified when the only quantity of interest being sought is the vector of approximate singular values. The resulting method is particularly effective for computing the nuclear norm of $A$, or more generally, other Schatten-$p$ norms. The report also describes how to compute an estimate of the errors incurred, at essentially negligible cost.

Equations

$$\text{rank }\mathsf{X}$$

$$\mathsf{X}\in \mathcal{C}$$

$$\begin{array}[]{ccccccccccc}\bm{\mathsf{A}}&=&\bm{\mathsf{U}}&\bm{\mathsf{T}}&\bm{\mathsf{V}}^{*}.\\ m\times n&&m\times m&m\times n&n\times n\end{array}$$

$$\Vert \mathsf{A}\Vert =\underset{\Vert \mathsf{x}\Vert =1}{\sup }\Vert \mathsf{A}\mathsf{x}{\Vert }_2={\sigma }_1(\mathsf{A}),\Vert \mathsf{A}{\Vert }_{\ast }=\sum _{i=1}^{\min (m,n)}{\sigma }_i(\mathsf{A}),\Vert \mathsf{A}{\Vert }_F={\left(\sum _{i=1}^m\sum _{j=1}^n{A}_{i,j}^2\right)}^{1/2}$$

$$\begin{tabular}[]{ccccccc}&&$b$&\ldelim[{2}{3mm}&$\bm{\mathsf{A}}_{11}$&$\bm{\mathsf{A}}_{12}$&\rdelim]{2}{3mm}\\ $\bm{\mathsf{A}}$&=&$m-b$&&$\bm{\mathsf{A}}_{21}$&$\bm{\mathsf{A}}_{22}$&\\ &&&&$b$&$n-b$&\end{tabular}.$$

$$\mathsf{T}={\mathsf{T}}_d+{\mathsf{T}}_u.$$

$$\sqrt{\sum_{i=1}^{\min(m,n)}\bigl{(}\sigma_{i}(\bm{\mathsf{A}})-\sigma_{i}(\bm{\mathsf{T}}_{d})\bigr{)}^{2}}\leq\|\bm{\mathsf{T}}_{u}\|_{F}.$$

$$\frac{|{\sigma }_i(\mathsf{A})-{\sigma }_i({\mathsf{T}}_d)|}{|{\sigma }_i(\mathsf{A})|},i=1,2,\dots ,\min (m,n).$$

Full text

Efficient nuclear norm approximation via the randomized UTV algorithm

N. Heavner111Department of Applied Mathematics, University of Colorado at Boulder, 526 UCB, Boulder, CO 80309-0526, USA, P.G. Martinsson111Department of Applied Mathematics, University of Colorado at Boulder, 526 UCB, Boulder, CO 80309-0526, USA

Abstract: The recently introduced algorithm randUTV provides a highly efficient technique for computing accurate approximations to all the singular values of a given matrix $\bm{\mathsf{A}}$. The original version of randUTV was designed to compute a full factorization of the matrix in the form $\bm{\mathsf{A}}=\bm{\mathsf{U}}\bm{\mathsf{T}}\bm{\mathsf{V}}^{*}$ where $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ are orthogonal matrices, and $\bm{\mathsf{T}}$ is upper triangular. The estimates to the singular values of $\bm{\mathsf{A}}$ appear along the diagonal of $\bm{\mathsf{T}}$. This manuscript describes how the randUTV algorithm can be modified when the only quantity of interest being sought is the vector of approximate singular values. The resulting method is particularly effective for computing the nuclear norm of $\bm{\mathsf{A}}$, or more generally, other Schatten-$p$ norms. The report also describes how to compute an estimate of the errors incurred, at essentially negligible cost.

1. Overview

This note describes an efficient algorithm for computing an accurate estimate for the nuclear norm $\|\cdot\|_{*}$ of a given matrix $\bm{\mathsf{A}}\in\mathbb{R}^{m\times n}$. The nuclear norm has recently found uses in numerical optimization, beginning with the introduction by Fazel et al.  in [3, 4] of its use as an effective heuristic for solving the rank minimization problem

[TABLE]

where $\bm{\mathsf{X}}\in\mathbb{R}^{m\times n}$ is the decision variable and $\mathcal{C}$ is some given convex constraint set. Recht et al. later proved in [11] that in certain cases, minimizing the nuclear norm also yields the theoretical solution to the corresponding rank minimization problem, solidifying the validity of the heuristic. The algorithm discussed in this note, randNN, may be used in a line search to choose the step size for nuclear norm minimization algorithms such as projected subgradient methods [11] or mirror descent [10, 2].

The recently proposed algorithm randUTV [8] is designed to compute, given a matrix $\bm{\mathsf{A}}\in\mathbb{R}^{m\times n}$, a factorization of the form

[TABLE]

In (1), the matrices $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ are orthogonal, $\bm{\mathsf{T}}$ is upper triangular, and the diagonal entries of $\bm{\mathsf{T}}$ are good approximations to the singular values of $\bm{\mathsf{A}}$. This note describes a modified version of randUTV which is of particular use in the sub-problem of conducting a line search during nuclear norm minimization. Specifically, we re-derive the randUTV algorithm from the perspective of nuclear norm estimation. In this case, since only the singular value estimates for $\bm{\mathsf{A}}$, rather than the entire matrix factorization, are desired, several steps of randUTV may be omitted or modified to yield randNN, which sees modest acceleration over randUTV and major acceleration over a full (SVD) computation of the singular values. We also mention a rough upper bound on the accuracy of the computed singular values that may computed as a part of randNN at very little extra cost from the middle matrix $\bm{\mathsf{T}}$.

The structure of this note is as follows. In Section 2, we review the notation used throughout the remainder. We derive the algorithm randNN in 3, and in 4 we discuss an error bound for the resulting approximations to the singular values. Finally, Section 5 contains numerical experiments exploring the accuracy of the singular values estimated by randNN and the performance of the algorithm compared to the industry standard.

2. Preliminaries

In this note, we use the notation $\bm{\mathsf{A}}\in\mathbb{R}^{m\times n}$ to denote a matrix $\bm{\mathsf{A}}$ of dimension $m\times n$ with real entries. $A_{i,j}$ denotes the $i$,$j$th entry of $\bm{\mathsf{A}}$. We use the notation of Golub and Van Loan [5] to specify submatrices: If $\bm{\mathsf{A}}$ is an $m\times n$ matrix, and $I=[i_{1},\,i_{2},\,\dots,\,i_{k}]$ and $J=[j_{1},\,j_{2},\,\dots,\,j_{\ell}]$ are index vectors, then $\bm{\mathsf{A}}(I,J)$ denotes the corresponding $k\times\ell$ submatrix. We let $\bm{\mathsf{A}}(I,:)$ denote the matrix $\bm{\mathsf{A}}(I,[1,\,2,\,\dots,\,n])$, and define $\bm{\mathsf{A}}(:,J)$ analogously. The transpose of a matrix $\bm{\mathsf{A}}$ is denoted $\bm{\mathsf{A}}^{*}$, and we use $\sigma_{i}(\bm{\mathsf{A}})$ to reference the $i$th leading singular value of $\bm{\mathsf{A}}$.

We measure matrices with either the spectral norm $\|\cdot\|,$ nuclear norm $\|\cdot\|_{*}$, or Frobenius norm $\|\cdot\|_{F}$, with definitions given by

[TABLE]

(where $\|\cdot\|_{2}$ is the Euclidean 2-norm for vectors in $\mathbb{R}^{n}$).

3. Description of the algorithm

3.1. Approach

Consider a matrix $\bm{\mathsf{A}}\in\mathbb{R}^{m\times n}$. Now partition $\bm{\mathsf{A}}$ as

[TABLE]

The construction of the algorithm will be guided by several fundamental observations:

• •

If $\bm{\mathsf{A}}_{12}=\bm{\mathsf{0}}$ and $\bm{\mathsf{A}}_{21}=\bm{\mathsf{0}}$, then $\|\bm{\mathsf{A}}\|_{*}=\|\bm{\mathsf{A}}_{11}\|_{*}+\|\bm{\mathsf{A}}_{22}\|_{*}$.

• •

Multiplying $\bm{\mathsf{A}}$ by an orthogonal matrix does not change the nuclear norm.

• •

Consider the matrix $\tilde{\bm{\mathsf{A}}}=\bm{\mathsf{U}}^{*}\bm{\mathsf{A}}\bm{\mathsf{V}}$, where $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ are orthogonal matrices whose first $b$ columns span the same space as the leading $b$ left and right singular vectors of $\bm{\mathsf{A}}$, respectively. Then giving $\tilde{\bm{\mathsf{A}}}$ the same partition as $\bm{\mathsf{A}}$, we have $\tilde{\bm{\mathsf{A}}}_{12}=\bm{\mathsf{0}}$ and $\tilde{\bm{\mathsf{A}}}_{21}=\bm{\mathsf{0}}$. For this reason, we will call such a $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ the “optimal” rotation matrices for the problem.

These remarks lay out a path for estimating $\|\bm{\mathsf{A}}\|_{*}$. We will choose $b$ to be a relatively small block size, say 50 or 100, and will compute an orthonormal matrices $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ whose first $b$ columns approximately span the same space as the respective leading $b$ left and right singular vectors of $\bm{\mathsf{A}}$. We will then form $\bm{\mathsf{A}}^{(1)}=\bm{\mathsf{U}}^{*}\bm{\mathsf{A}}\bm{\mathsf{V}}$, after which we will have $\bm{\mathsf{A}}^{(1)}_{21}=\bm{\mathsf{0}}$ and $\|\bm{\mathsf{A}}^{(1)}_{12}\|$ is small. At this point, we may estimate the norm by finding the norms of the small $\bm{\mathsf{A}}^{(1)}_{11}$ block and the large $\bm{\mathsf{A}}^{(1)}_{22}$ block separately. $\|\bm{\mathsf{A}}^{(1)}_{11}\|_{*}$ may be computed efficiently with an SVD computation, and we may recursively apply our strategy for $\bm{\mathsf{A}}$ to $\bm{\mathsf{A}}^{(1)}_{22}$ to estimate its norm.

3.2. Computing $\bm{\mathsf{V}}$

At this point in the derivation, the most pressing question is how to find $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ that satisfy the aforementioned conditions. We first consider the construction of $\bm{\mathsf{V}}$, at which point $\bm{\mathsf{U}}$ will easily follow. Recent work in randomized subspace iteration, discussed in papers including [7, 12, 14, 6], enables the efficient computation of the desired $\bm{\mathsf{V}}$.

In particular, consider a sampling matrix $\bm{\mathsf{Y}}=\bm{\mathsf{A}}^{*}\bm{\mathsf{G}}$, where $\bm{\mathsf{G}}$ is an $m\times b$ Gaussian random matrix. Then with high probability, the column space of $\bm{\mathsf{Y}}$ often aligns closely to the subspace spanned by the leading $b$ right singular vectors of $\bm{\mathsf{A}}$. As discussed in [12, 7], in certain situations when the decay of singular values of $\bm{\mathsf{A}}$ makes the alignment suboptimal, then using $\bm{\mathsf{Y}}=(\bm{\mathsf{A}}^{*}\bm{\mathsf{A}})^{q}\bm{\mathsf{A}}^{*}\bm{\mathsf{G}}$ instead provides considerable correction, where $q$ is some small nonnegative integer. Typically, $q=0,1,2$ suffices to obtain a close-to-optimal approximation to the desired subspace.

Thus, the following steps compute a matrix $\bm{\mathsf{V}}$ whose columns span approximately the same subspace as the leading $b$ right singular vectors of $\bm{\mathsf{A}}$:

• •

Draw an $n\times b$ random Gaussian matrix $\bm{\mathsf{G}}$.

• •

Compute $\bm{\mathsf{Y}}=(\bm{\mathsf{A}}^{*}\bm{\mathsf{A}})^{q}\bm{\mathsf{A}}^{*}\bm{\mathsf{G}}$ for some small nonnegative integer $q$.

• •

Compute an orthonormal basis for the column space of $\bm{\mathsf{Y}}$ with a QR factorization using Householder reflectors to obtain $\bm{\mathsf{Y}}=\bm{\mathsf{V}}\bm{\mathsf{R}}$.

3.3. Computing $\bm{\mathsf{U}}$

To compute $\bm{\mathsf{U}}$, we first note that if our computed $\bm{\mathsf{V}}$ were optimal, then the first $b$ columns of $\bm{\mathsf{A}}\bm{\mathsf{V}}$ would be composed exclusively of linear combinations of the leading $b$ left singular vectors of $\bm{\mathsf{A}}$. Thus, to form $\bm{\mathsf{U}}$, we may simply compute an orthonormal basis for the first $b$ columns of $\bm{\mathsf{A}}\bm{\mathsf{V}}$ via Householder reflectors to obtain $\bm{\mathsf{A}}\bm{\mathsf{V}}(:,1:b)=\bm{\mathsf{U}}\bm{\mathsf{R}}$.

Remark 1**.**

In applications, the input matrix $\bm{\mathsf{A}}$ is often rank-deficient, with numerical rank $k$ determined by the $k$ for which $\sigma_{i}(\bm{\mathsf{A}})$ is less than some specified tolerance for $k<i\leq n$. In the case where $k\ll n$, randNN may be sped up substantially by terminating the algorithm once the estimated singular values drop below some user-defined threshold.

4. Error bounds

In this section, we make note of an error bound for the accuracy of the estimated singular values resulting from the randNN algorithm.

First, note that after applying $\lceil n/b\rceil$ left and right rotation matrices to $\bm{\mathsf{A}}$ as described in 3, we are left with an upper triangular matrix which we shall call $\bm{\mathsf{T}}$ whose mass is concentrated in blocks of size $b\times b$ along its diagonal. Therefore, let $\bm{\mathsf{T}}_{d}$ be the block diagonal matrix consisting of the $b\times b$ main diagonal blocks of $\bm{\mathsf{T}}$, and let $\bm{\mathsf{T}}_{u}$ be the upper triangular matrix defined by the relation

[TABLE]

Observe that the estimated singular values of $\bm{\mathsf{A}}$ in our algorithm consist precisely of the computed singular values of $\bm{\mathsf{T}}_{d}$. Next, since $\bm{\mathsf{A}}$ and $\bm{\mathsf{T}}$ only differ by orthogonal transforms, we have that $\|\bm{\mathsf{A}}\|_{*}=\|\bm{\mathsf{T}}\|_{*}$, and more specifically, $\sigma_{i}(\bm{\mathsf{A}})=\sigma_{i}(\bm{\mathsf{T}})$ for $i=1,2,\ldots,\min(m,n)$. Thus an error bound given in [13] and [9] states that

[TABLE]

Since the calculation of $\|\bm{\mathsf{T}}_{u}\|_{F}$ is $\mathcal{O}(n^{2})$, this bound may be calculated at negligible relative cost as part of randNN, serving as an assurance of the estimation’s validity.

5. Numerical experiments

5.1. Computational speed

In this section, we investigate the speed of the proposed algorithm randNN and compare it to the speed of an exact computation of the singular values using LAPACK’s highly optimized dgesvd.

All experiments reported in this note were performed on an Intel Core i7-6700K processor (4.0 GHz) with 4 cores. In order to be able to show scalability results, the clock speed was throttled at 4.0 GHz, turning off so-called turbo boost. Other details of interest include that the OS used was Ubuntu (Version 16.04.2), and the code was compiled with Intel’s icc (Version 17.0.4.196). Main routines for computing the singular values (dgesvd) were taken from Intel’s MKL library (Version 2017.0.3). Standard BLAS routines used in our implementation of randNN were also taken from the Intel MKL library.

Each of the three algorithms we tested (randNN,randUTV,SVD) was applied to double-precision real matrices of size $n\times n$. We report the following times:

• •

$T_{\text{svd}}$. The time in seconds for the LAPACK function dgesvd from Intel’s MKL, where the orthogonal matrices $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ in the SVD were not built.

• •

$T_{\text{randUTV}}$. The time in seconds for the function randUTV described in [8]. The authors’ original implementation was used, and the orthogonal matrices $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ in the UTV decomposition were not built.

• •

$T_{\text{randNN}}$. The time in seconds for our implementation of randNN.

In all cases, we used a block size of $b=64$. While likely not optimal for all problem sizes, this block size yields near best performance and allows us to easily compare and contrast the performance of the different implementations.

As expected, Figure 2 shows that both randUTV and randNN are markedly faster than LAPACK’s dgesvd, even in the single case. When more cores are added to the computation, the gain in speed is even more dramatic, hi-lighting the fact that the SVD algorithm is not designed to make efficient use of parallel computing architectures. Finally the increase in speed from randUTV to randNN is modest, but in a line search setting where the nuclear norm must be computed many times, this saving can add up over the course of the solution of the rank minimization problem.

5.2. Errors

In this section, we report the results of numerical experiments that were conducted to test the accuracy of the approximation to the singular values provided by randNN. Specifically, we compare the estimated singular values to the true singular values of two different test matrices:

• •

Matrix 1 (S shaped decay): This is an $n\times n$ matrix of the form $\bm{\mathsf{A}}=\bm{\mathsf{U}}\bm{\mathsf{D}}\bm{\mathsf{V}}^{*}$ where $\bm{\mathsf{U}}$ and $\bm{\mathsf{V}}$ are randomly drawn matrices with orthonormal columns (obtained by performing QR on a random Gaussian matrix), and where the diagonal entries of $\bm{\mathsf{D}}$ are chosen to first hover around 1, then decay rapidly, and then level out at $10^{-6}$, as shown in Figure 3 (black line) on the left.

• •

Matrix 2 (Single Layer BIE): This matrix is the result of discretizing a Boundary Integral Equation (BIE) defined on a smooth closed curve in the plane. To be precise, we discretized the so called “single layer” operator associated with the Laplace equation using a $6^{\rm th}$ order quadrature rule designed by Alpert [1]. This is a well-known ill-conditioned problem for which column pivoting is essential in order to stably solve the corresponding linear system.

In each case, a matrix size of $m=n=5000$ was used with a “power iteration parameter” (see Section 3.2) of $q=2$ and a block size of $b=64$.

For each test matrix, we plot both the estimated and true singular values themselves and the relative error of the estimates

[TABLE]

6. Availability of code

An implementation of the discussed algorithm is available under 3-clause (modified) BSD license from:

[TABLE]