Efficient nuclear norm approximation via the randomized UTV algorithm
Nathan Heavner Per-Gunnar Martinsson












widequote
      
      
  










thmTheorem
lemma[thm]Lemma
prop[thm]Proposition
corollary[thm]Corollary
conj[thm]Conjecture
definition
remarkRemark
exerciseExercise
definitionDefinition































Efficient nuclear norm approximation via the randomized UTV algorithm




N. Heavner[1],
P.G. Martinsson[1]





140mm
Abstract: The recently introduced algorithm  provides a highly efficient technique for computing accurate approximations to all the singular values of a given matrix A. The original version of  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  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.





[1]Department of Applied Mathematics, University of Colorado at Boulder, 526 UCB, Boulder, CO 80309-0526, USA








§ OVERVIEW



This note describes an efficient algorithm for computing an accurate estimate for the nuclear norm ·_* of a given matrix A∈ℝ^m × n. The nuclear norm has recently found uses in numerical optimization, beginning with the introduction by Fazel et al.  in <cit.> of its use as an effective heuristic for solving the rank minimization problem

    minimize   =rank X
    subject to   =X∈𝒞

where X∈ℝ^m × n is the decision variable and 𝒞 is some given convex constraint set. Recht et al. later proved in <cit.> 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, , may be used in a line search to choose the step size for nuclear norm minimization algorithms such as projected subgradient methods <cit.> or mirror descent <cit.>.

The recently proposed algorithm  <cit.> is designed to compute, given a  matrix A∈ℝ^m × n, a factorization of the form

    [    A    =    U    T V^*.; m× n m× m m× n n× n ]

In (<ref>), the matrices U and V are orthogonal, T is upper triangular, and the diagonal entries of T are good approximations to the singular values of A. This note describes a modified version of  which is of particular use in the sub-problem of conducting a line search during nuclear norm minimization. Specifically, we re-derive the  algorithm from the perspective of nuclear norm estimation. In this case, since only the singular value estimates for A, rather than the entire matrix factorization, are desired, several steps of  may be omitted or modified to yield , which sees modest acceleration over  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  at very little extra cost from the middle matrix T.

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




§ PRELIMINARIES




In this note, we use the notation A∈ℝ^m × n to denote a matrix A of dimension m × n with real entries. A_i,j denotes the i,jth entry of A. We use the notation of Golub and Van Loan <cit.> to specify submatrices: If A is an m× n matrix, and I = [i_1, i_2, …, i_k] and J = [j_1, j_2, …, j_ℓ] are index vectors, then A(I,J) denotes the corresponding k ×ℓ submatrix. We let A(I,:) denote the matrix A(I,[1, 2, …, n]), and define A(:,J) analogously. The transpose of a matrix A is denoted A^*, and we use σ_i(A) to reference the ith leading singular value of A.

We measure matrices with either the spectral norm ·, nuclear norm ·_*, or Frobenius norm ·_F, with definitions given by

    A = sup_x⃗ = 1Ax⃗_2 = σ_1(A),   A_* =  ∑_i=1^min(m,n)σ_i(A),   A_F =  (  ∑_i=1^m ∑_j=1^n A_i,j^2   )^1/2

(where ·_2 is the Euclidean 2-norm for vectors in ℝ^n).




§ DESCRIPTION OF THE ALGORITHM







 §.§ Approach



Consider a matrix A∈ℝ^m × n. Now partition A as


    b   [23mm   A_11   A_12   ]23mm
    A    =    m - b      A_21   A_22   
                b   n - b   .


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



  * If A_12 = 0 and A_21 = 0, then A_* = A_11_* + A_22_*.

  * Multiplying A by an orthogonal matrix does not change the nuclear norm.

  * Consider the matrix Ã = U^* AV, where U and V are orthogonal matrices whose first b columns span the same space as the leading b left and right singular vectors of A, respectively. Then giving Ã the same partition as A, we have Ã_12 = 0 and Ã_21 = 0. For this reason, we will call such a U and V  the “optimal” rotation matrices for the problem.


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





 §.§ Computing V




At this point in the derivation, the most pressing question is how to find U and V that satisfy the aforementioned conditions. We first consider the construction of V, at which point U will easily follow. Recent work in randomized subspace iteration, discussed in papers including <cit.>, enables the efficient computation of the desired V.

In particular, consider a sampling matrix Y = A^* G, where G is an m × b Gaussian random matrix. Then with high probability, the column space of  Y often aligns closely to the subspace spanned by the leading b right singular vectors of A. As discussed in <cit.>, in certain situations when the decay of singular values of A makes the alignment suboptimal, then using Y = (A^* A)^q A^* 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 V whose columns span approximately the same subspace as the leading b right singular vectors of A:



  * Draw an n × b random Gaussian matrix G.

  * Compute Y  = (A^* A)^q A^* G for some small nonnegative integer q.

  * Compute an orthonormal basis for the column space of Y with a QR factorization using Householder reflectors to obtain Y = VR.





 §.§ Computing U



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


In applications, the input matrix A is often rank-deficient, with numerical rank k determined by the k for which σ_i(A) is less than some specified tolerance for k < i ≤ n. In the case where k ≪ n,  may be sped up substantially by terminating the algorithm once the estimated singular values drop below some user-defined threshold.








§ ERROR BOUNDS




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

First, note that after applying ⌈ n/b ⌉ left and right rotation matrices to A as described in <ref>, we are left with an upper triangular matrix which we shall call T whose mass is concentrated in blocks of size b × b along its diagonal. Therefore, let T_d be the block diagonal matrix consisting of the b × b main diagonal blocks of T, and let T_u be the upper triangular matrix defined by the relation

    T = T_d + T_u.

Observe that the estimated singular values of A in our algorithm consist precisely of the computed singular values of T_d. Next, since A and T only differ by orthogonal transforms, we have that A_* = T_*, and more specifically, σ_i(A) = σ_i(T) for i=1,2,…,min(m,n). Thus an error bound given in <cit.> and <cit.> states that

    √(∑_i=1^min(m,n)(σ_i(A) - σ_i(T_d))^2)≤T_u _F.


Since the calculation of T_u _F is 𝒪(n^2), this bound may be calculated at negligible relative cost as part of , serving as an assurance of the estimation's validity.




§ NUMERICAL EXPERIMENTS







 §.§ Computational speed



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

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  (Version 17.0.4.196). Main routines for computing the singular values () were taken from Intel's MKL library (Version 2017.0.3). Standard BLAS routines used in our implementation of  were also taken from the Intel MKL library.

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


  * T_svd. The time in seconds for the LAPACK function  from Intel's MKL, where the orthogonal matrices U and V in the SVD were not built.

  * T_randUTV. The time in seconds for the function  described in <cit.>. The authors' original implementation was used, and the orthogonal matrices U and V in the UTV decomposition were not built.

  * T_randNN. The time in seconds for our implementation of .


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 <ref> shows that both  and  are markedly faster than LAPACK's , 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  to  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.




 §.§ 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 . 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 × n matrix of the form A = UDV^* where U and V are randomly drawn matrices with orthonormal columns (obtained by performing QR on a random Gaussian matrix), and where the diagonal entries of D are chosen to first hover around 1, then decay rapidly, and then level out at 10^-6, as shown in Figure <ref> (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^ th order quadrature rule designed by Alpert <cit.>. 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 <ref>) 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

    | σ_i(A) - σ_i(T_d) |/| σ_i(A) |,    i=1,2,…,min(m,n).







§ AVAILABILITY OF CODE



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

    https://github.com/nheavner/nn_code



amsplain




