# splitting numerical integration for matrix completion

Qianqian Song

###### Abstract

Low rank matrix approximation is a popular topic in machine learning. In this
paper, we propose a new algorithm for this topic by minimizing the least-
squares estimation over the Riemannian manifold of fixed-rank matrices. The
algorithm is an adaptation of classical gradient descent within the framework
of optimization on manifolds. In particular, we reformulate an unconstrained
optimization problem on a low-rank manifold into a differential dynamic
system. We develop a splitting numerical integration method by applying a
splitting integration scheme to the dynamic system. We conduct the convergence
analysis of our splitting numerical integration algorithm. It can be
guaranteed that the error between the recovered matrix and true result is
monotonically decreasing in the Frobenius norm. Moreover, our splitting
numerical integration can be adapted into matrix completion scenarios.
Experimental results show that our approach has good scalability for large-
scale problems with satisfactory accuracy.

##  1 Introduction

Computing an efficient and reliable low-rank approximation of a given matrix
is a fundamental task in many machine learning problems, such as principal
component analysis Jolliffe (2002), face recognition Muller et al. (2004) and
large scale data compression Drineas et al. (2006); Huang and others (2020).
It is well-known that the truncated singular value decomposition (SVD)
provides the best low-rank approximation to the matrix in question.

Specifically, for a given matrix $\mathbf{M}\in\mathbb{R}^{m\times n}$ with
$m\geq n$, its SVD is defined as

$\mathbf{M}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T},$

where $\mathbf{U}=[\mathbf{u}_{1},\ldots,\mathbf{u}_{n}]$ is an $m\times n$
column orthonormal matrix, $\mathbf{V}=[\mathbf{v}_{1},\ldots,\mathbf{v}_{n}]$
is an orthonormal matrix, and
$\mathbf{\Sigma}=\mathrm{diag}(\sigma_{1},\ldots,\sigma_{n})$ is a diagonal
matrix with the diagonal entries
$\sigma_{1}\geq\sigma_{2}\geq\cdots\sigma_{n}\geq 0$. Moreover,
$\mathbf{u}_{j}$ and $\mathbf{v}_{j}$ are called the left and right singular
vectors corresponding to $\sigma_{j}$, the $j$-th largest singular value of
$\mathbf{M}$. For any $1\leq r\leq n$, then

$\mathbf{M}_{r}=[\mathbf{u}_{1},\ldots,\mathbf{u}_{r}]\mathrm{diag}(\sigma_{1},\ldots,\sigma_{r})[\mathbf{v}_{1},\ldots,\mathbf{v}_{r}]^{T}$

is the truncated SVD of $\mathbf{M}$ of rank at most $k$, which is unique only
if $\sigma_{r+1}<\sigma_{r}$. The assumption that $m\geq n\gg r$ will be
maintained throughout this paper for clarity of statement. All results will
also hold for $m<n$, applied on $\mathbf{M}^{T}$.

It is well known that $\mathbf{M}_{r}$ is the best rank-$r$ approximation to
$\mathbf{M}$ Eckart and Young (1936); Golub and Van Loan (2012). The truncated
SVD can be cast into a fixed-rank optimzation problem. That is,

$\displaystyle\underset{\mathbf{Y}}{\arg\min}\|\mathbf{M}-\mathbf{Y}\|_{F},$
(1) $\displaystyle\text{s.t.}\ \text{rank}(\mathbf{Y})=r.$

A bunch of algorithms, such as the Lanczos algorithm Golub and Van Loan
(2012), randomized SVD algorithm Halko et al. (2011), and subspace iteration
Gu (2015), have been proposed to solve this problem.

An equivalent formulation is given by introducing the concept of matrix
manifold as follows:

$\displaystyle\underset{\mathbf{Y}}{\arg\min}\frac{1}{2}\|\mathbf{M}-\mathbf{Y}\|_{F}^{2},$
(2) $\displaystyle\text{s.t.}\ \mathbf{Y}\in\mathcal{M}_{r},$

where $\mathcal{M}_{r}$ represents a rank-$r$ smooth manifold. This
formulation can be viewed as an unconstrained optimization problem on the low-
rank manifold. Accordingly, Koch and Lubich Koch and Lubich (2007) proposed a
dynamical low rank approximation approach, which is a gradient descent
procedure on the manifold in essence.

In this paper, we attempt to make the gradient descent process “finer-
grained.” In particular, we view the problem from a perspective of
differential dynamic systems and employ a splitting scheme to solve it.
Accordingly, we devise a novel method for low rank matrix approximation that
we call a splitting numerical integration method. Theoretical analysis
guarantees that the splitting numerical integration algorithm converges
asymptotically. In addition, we apply splitting numerical integration to
matrix completion scenarios, where the matrix $\mathbf{M}$ is partially
observed. Empirical results are also encouraging, especially on large-scale
datasets.

The remainder of the paper is organized as follows. Section 2 presents the
notation frequently used in this paper and the problem formulation. Section 3
describes our algorithm and theoretical analysis. Empirical results are given
in Section 4.

##  2 Notation and Preliminaries

First of all, we present the notation used in this paper. Let $\mathbf{I}_{m}$
be the $m\times m$ identity matrix. Given a matrix
$\mathbf{Y}\in\mathbb{R}^{m\times n}$, $\|\mathbf{Y}\|_{F}$ denotes the
Frobenius norm of $\mathbf{Y}$ and $\|\mathbf{Y}\|_{2}$ denotes the spectral
norm.

It is well established that every rank-$r$ matrix
$\mathbf{Y}\in\mathbb{R}^{m\times n}$ can be written in the form

$\displaystyle\mathbf{Y}=\mathbf{U}\mathbf{S}\mathbf{V}^{T},$ (3)

where $\mathbf{U}\in\mathbb{R}^{m\times r}$ and
$\mathbf{V}\in\mathbb{R}^{n\times r}$ are column orthonormal, i.e.,
$\mathbf{U}^{T}\mathbf{U}=\mathbf{I}_{r}\;\mbox{ and
}\mathbf{V}^{T}\mathbf{V}=\mathbf{I}_{r}$, and
$\mathbf{S}\in\mathbb{R}^{r\times r}$ is nonsingular. Notice that here we do
not require $\mathbf{S}$ to be the diagonal matrix of the singular values. The
representation in Eqn. (3) is not unique because
$\mathbf{Y}=\hat{\mathbf{U}}\hat{\mathbf{S}}\hat{\mathbf{V}}^{T}$ is another
representation where $\hat{\mathbf{U}}=\mathbf{U}\mathbf{P}$,
$\hat{\mathbf{V}}=\mathbf{V}\mathbf{Q}$, and
$\hat{\mathbf{S}}=\mathbf{P}^{T}\mathbf{S}\mathbf{Q}$ whenever
$\mathbf{P},\mathbf{Q}\in\mathbb{R}^{r\times r}$ are any orthonormal matrices.

As a substitute for the non-uniqueness in Eqn. (3), we will use a unique
decomposition in the tangent space. Let $\mathcal{V}_{m,r}$ represent the
Stiefel manifold of real column orthonormal matrices of size $m\times r$
($m>r$). The tangent space at the point $\mathbf{U}\in\mathcal{V}_{m,r}$ is
defined as:

$\displaystyle\mathcal{T}_{\mathbf{U}}\mathcal{V}_{m,r}$
$\displaystyle=\\{\delta\mathbf{U}\in\mathbb{R}^{m\times
r}:\delta\mathbf{U}^{T}\mathbf{U}+\mathbf{U}^{T}\delta\mathbf{U}=\mathbf{0}\\}$
$\displaystyle=\\{\delta\mathbf{U}\in\mathbb{R}^{m\times
r}:\mathbf{U}^{T}\delta\mathbf{U}\in so(r)\\},$

where $so(r)$ denotes the space of skew-symmetric real $r\times r$ matrices.
Consider the extended tangent map of
$(\mathbf{S},\mathbf{U},\mathbf{V})\longmapsto\mathbf{Y}=\mathbf{U}\mathbf{S}\mathbf{V}^{T}$,

$\displaystyle\mathbb{R}^{r\times
r}\times\mathcal{T}_{\mathbf{U}}\mathcal{V}_{m,r}\times\mathcal{T}_{\mathbf{V}}\mathcal{V}_{n,r}$
$\displaystyle\xrightarrow{}\mathcal{T}_{\mathbf{Y}}M_{r}\times so(r)\times
so(r),$ $\displaystyle(\delta\mathbf{S},\delta\mathbf{U},\delta\mathbf{V})$
$\displaystyle\xrightarrow{}(\delta\mathbf{U}\mathbf{S}\mathbf{V}^{T}+\mathbf{U}\delta\mathbf{S}\mathbf{V}^{T}+\mathbf{U}\mathbf{S}\delta\mathbf{V}^{T},\mathbf{U}^{T}\delta\mathbf{U},\mathbf{V}^{T}\delta\mathbf{V}).$

The manifold of rank-$r$ matrices is denoted by $\mathcal{M}_{r}$. The tangent
space at any $\mathbf{Y}\in\mathcal{M}_{r}$ is denoted by
$\mathcal{T}_{\mathbf{Y}}\mathcal{M}_{r}$, which is defined as follows. Every
$\delta\mathbf{Y}\in\mathcal{T}_{\mathbf{Y}}\mathcal{M}_{r}$ can be written
into the following form:

$\displaystyle\delta{\mathbf{Y}}=\delta{\mathbf{U}}\mathbf{S}\mathbf{V}^{T}+\mathbf{U}\delta{\mathbf{S}}\mathbf{V}^{T}+\mathbf{U}\mathbf{S}\delta{\mathbf{V}}^{T},$
(4)

where $\delta\mathbf{S}\in\mathbb{R}^{r\times r}$,
$\delta\mathbf{U}\in\mathcal{T}_{\mathbf{U}}\mathcal{V}_{m,r}$ and
$\delta\mathbf{V}\in\mathcal{T}_{\mathbf{V}}\mathcal{V}_{n,r}$. Furthermore,
$\delta\mathbf{S}$, $\delta\mathbf{U}$ and $\delta\mathbf{V}$ are uniquely
determined by $\delta\mathbf{Y}$ if we impose the orthogonality constraints:

$\displaystyle\mathbf{U}^{T}\delta\mathbf{U}=\mathbf{0},$ (5)
$\displaystyle\mathbf{V}^{T}\delta\mathbf{V}=\mathbf{0}.$

The projection operators onto the spaces spanned by the columns of
$\mathbf{U}$ and $\mathbf{V}$, and their orthogonal complements are defined as

$\displaystyle P_{\mathbf{U}}=\mathbf{U}\mathbf{U}^{T},$ $\displaystyle
P_{\mathbf{V}}=\mathbf{V}\mathbf{V}^{T},$ $\displaystyle
P^{\bot}_{\mathbf{U}}=\mathbf{I}_{m}-\mathbf{U}\mathbf{U}^{T},$ $\displaystyle
P^{\bot}_{\mathbf{V}}=\mathbf{I}_{n}-\mathbf{V}\mathbf{V}^{T}.$

Then $\delta\mathbf{S}$, $\delta\mathbf{U}$ and $\delta\mathbf{V}$ are
uniquely determined by $\delta\mathbf{Y}$ as follows:

$\displaystyle\delta\mathbf{S}=\mathbf{U}^{T}\delta\mathbf{Y}\mathbf{V},$
$\displaystyle\delta\mathbf{U}=P^{\bot}_{\mathbf{U}}\delta\mathbf{Y}\mathbf{V}\mathbf{S}^{-1},$
$\displaystyle\delta\mathbf{V}=P^{\bot}_{\mathbf{V}}\delta\mathbf{Y}^{T}\mathbf{U}\mathbf{S}^{-T}.$

###  2.1 Dynamic Low Rank Approximation

Let us return to Problem (2) and let

$\displaystyle f(\mathbf{Y})=\frac{1}{2}\|\mathbf{Y}-\mathbf{M}\|_{F}^{2}$

denote the objective function. The gradient of $f(\mathbf{Y})$ in $\mathbf{Y}$
can be written as:

$\displaystyle\nabla f(\mathbf{Y})=\mathbf{Y}-\mathbf{M}.$

Recall that the constraint $\mathbf{Y}\in\mathcal{M}_{r}$, Riemannian
gradient, denoted $\nabla\mathbf{Y}$, is used instead of $\nabla
f(\mathbf{Y})$, which is a specific tangent vector corresponding to the
direction of steepest ascent of $f(\mathbf{Y})$ but restricted to the tangent
space $\mathcal{T}_{\mathbf{Y}}\mathcal{M}_{r}$. It can be solved via the
following optimization problem:

$\displaystyle\underset{\nabla{\mathbf{Y}}\in\mathbb{R}^{m\times
n}}{\arg\min}\;\frac{1}{2}\|\nabla{\mathbf{Y}}-\nabla{\mathbf{A}}\|_{F}^{2}$
(6) $\displaystyle\text{s.t.}\
\nabla{\mathbf{Y}}\in\mathcal{T}_{X}\mathcal{M}_{r},$

where we denote $\nabla{\mathbf{A}}\triangleq-\nabla f(\mathbf{Y})$ for
notational simplicity and $\nabla\mathbf{A}$ can be seen as a given constant
at each iteration. The solution of Problem (6) is well-studied in the
differential geometry literature. That is,

######  Proposition 1.

_Koch and Lubich (2007)_ Let
$\mathbf{Y}=\mathbf{U}\mathbf{S}\mathbf{V}^{T}\in\mathcal{M}_{r}$ where
$\mathbf{S}\in\mathbb{R}^{r\times r}$ is nonsingular, and
$\mathbf{U}\in\mathbb{R}^{m\times r}$ and $\mathbf{V}\in\mathbb{R}^{n\times
r}$ are column orthonormal. Then the solution to Problem (6) can be written in
the following form:

$\displaystyle\nabla{\mathbf{Y}}=\nabla{\mathbf{U}}\mathbf{S}\mathbf{V}^{T}+\mathbf{U}\nabla{\mathbf{S}}\mathbf{V}^{T}+\mathbf{U}\mathbf{S}\nabla{\mathbf{V}}^{T},$
(7)

where

$\displaystyle\nabla{\mathbf{S}}=\mathbf{U}^{T}\nabla{\mathbf{A}}\mathbf{V},$
(8)
$\displaystyle\nabla{\mathbf{U}}=P_{\mathbf{U}}^{\bot}\nabla{\mathbf{A}}\mathbf{V}\mathbf{S}^{-1},$
$\displaystyle\nabla{\mathbf{V}}=P_{\mathbf{V}}^{\bot}\nabla{\mathbf{A}}^{T}\mathbf{U}\mathbf{S}^{-T}.$

The resulting algorithm is a simple gradient descent procedure. That is,

$\displaystyle\mathbf{U}\Leftarrow\mathbf{U}+\epsilon\nabla\mathbf{U},\;\mathbf{S}\Leftarrow\mathbf{S}+\epsilon\nabla\mathbf{S},\;\mbox{
and }\;\mathbf{V}\Leftarrow\mathbf{V}+\epsilon\nabla\mathbf{V},$ (9)

where $\epsilon$ is stepsize. The algorithm was firstly proposed in Koch and
Lubich (2007) and called dynamical low-rank approximation, a landmark in
solving SVD in the view of dynamic system. Prior workof solving SVD in dynamic
system is about full SVD, without considering truncated SVD.

According to the low rank assumption, the matrix $\mathbf{Y}$ does not come up
explicitly for ease of computation. Instead, $\mathbf{Y}$ is represented via
the product of $\mathbf{U},\mathbf{S}$ and $\mathbf{V}^{T}$ (at order of
$O(mr+nr)$), as shown in Eqn. (3). Similarly, the Riemannian gradient
$\nabla\mathbf{Y}$ can be also represented by using relatively small matrices
(at order of $O(mr+nr)$) as shown in Eqn. (7). However, for notational
convenience, we use $\mathbf{Y}$ and $\nabla\mathbf{Y}$ instead of
$\mathbf{U}\mathbf{S}\mathbf{V}^{T}$ and
$(\nabla\mathbf{U})\mathbf{S}\mathbf{V}^{T}+\mathbf{U}(\nabla\mathbf{S})\mathbf{V}^{T}+\mathbf{U}\mathbf{S}(\nabla\mathbf{V})^{T}$,
respectively.

Though the extant dynamical low rank approximation is a benchmark work in
solving low rank approximation in view of dynamic system, it is not
competitive compared with state-of-the-art low rank approximation approaches
due to the following issues.

First, this framework fails to exploit geometric information sufficiently.
Specifically, on the fixed-rank manifold, from
$\mathbf{Y}_{i}=\mathbf{U}_{i}\mathbf{S}_{i}\mathbf{V}_{i}^{T}$ to
$\mathbf{Y}_{i+1}=\mathbf{U}_{i+1}\mathbf{S}_{i+1}\mathbf{V}_{i+1}^{T}$, the
algorithm only makes use of the Riemannian gradient of $\mathbf{Y}_{i}$.
Second, since the discretization error of the dynamic system is proportional
to the stepsize Koch and Lubich (2007), the stepsize is required to be
extremely small, which limits the convergence rate. Third, some extra
operations are introduced to ensure the column orthogonality of $\mathbf{U}$
and of $\mathbf{V}$, which is deemed to be a brutal strategy, often leading to
loss of information Koch and Lubich (2007).

To address these issues, we seek a widely used scheme in differential systems
to compute gradient descent. Specifically, we use a splitting integration
technique to update three components ($\mathbf{U}$, $\mathbf{S}$,
$\mathbf{V}$) step-by-step, making a more sufficient use of geometric
information.

##  3 Methodology

In this section we present our method. We first give a novel view of dynamic
system for Problem (6). With this view, we use a splitting integration scheme
to devise a dynamic flow subspace method for solving Problem (2). Finally, we
give the convergence analysis of the method and extend it into full SVD and
low rank matrix completion.

###  3.1 Dynamic System

Departing from a perspective of differential systems, we study the solution of
the optimization problem in (6). In particular, we consider an alternative
formulation for Problem (2) in the form of a dynamic system:

$\displaystyle\underset{{\dot{\mathbf{Y}}\in\mathbb{R}^{m\times
n}}}{\arg\min}\;\frac{1}{2}\|\dot{\mathbf{Y}}-\dot{\mathbf{A}}\|_{F}^{2}$ (10)
$\displaystyle\text{s.t.}\ \dot{\mathbf{Y}}\in T_{X}\mathcal{M}_{r},$

where $\mathbf{Y}$ is regarded as a time-dependent matrix such that
$\mathbf{Y}=\mathbf{Y}(t)$ and $\dot{\mathbf{Y}}$ denotes the derivative of
$\mathbf{Y}$ w.r.t. time, and
$\dot{\mathbf{A}}\triangleq\nabla\mathbf{A}=-\nabla f(\mathbf{Y})$ in our
case. Notice that it is the continuous version of Problem (6).

According to the Galerkin condition on the tangent space
$\mathcal{T}_{\mathbf{Y}}\mathcal{M}_{r}$ in numerical analysis Hairer et al.
(2006), Problem (10) is equivalent to the following projection:

$\displaystyle\text{finding}\
\dot{\mathbf{Y}}\in\mathcal{T}_{\mathbf{Y}}\mathcal{M}_{r}\ \text{such that}$
(11)
$\displaystyle\left\langle\dot{\mathbf{Y}}-\dot{\mathbf{A}},\delta\mathbf{Y}\right\rangle=0\
\ \text{for\ all\
}\delta\mathbf{Y}\in\mathcal{T}_{\mathbf{Y}}\mathcal{M}_{r}.$

Furthermore, Problem (11) can be transformed into the following form:

$\displaystyle\dot{\mathbf{Y}}=\tilde{P}_{\mathbf{Y}}(\dot{\mathbf{A}}),$ (12)

where $\tilde{P}_{\mathbf{Y}}(\cdot)$ is a projection operator, defined as

$\displaystyle\tilde{P}_{\mathbf{Y}}(\mathbf{B})=\mathbf{B}-\tilde{P}_{\mathbf{Y}}^{\bot}(\mathbf{B})$
(13) $\displaystyle\text{with}\ \
\tilde{P}_{\mathbf{Y}}^{\bot}(\mathbf{B})=P^{\bot}_{\mathbf{U}}\mathbf{B}P^{\bot}_{\mathbf{V}}\
\;\ \forall\ \mathbf{B}\in\mathbb{R}^{m\times n}.$

Substituting Eqn. (13) into Eqn. (12), we have the following dynamic system:

$\displaystyle\dot{\mathbf{Y}}=\mathbf{U}\mathbf{U}^{T}\dot{\mathbf{A}}-\mathbf{U}\mathbf{U}^{T}\dot{\mathbf{A}}\mathbf{V}\mathbf{V}^{T}+\dot{\mathbf{A}}\mathbf{V}\mathbf{V}^{T}.$
(14)

###  3.2 Dynamic Flow Subspace Method

Our current concern is to solve the differential equation in (14). We resort
to a splitting scheme Hairer et al. (2006). In particular, let $\mathcal{L}$
be a local generator Leimkuhler and Matthews (2013) corresponding to the exact
solution to Problem (14) and separate it into several sub-generators as
follows:

$\displaystyle\mathcal{L}=\mathcal{L}_{A}+\mathcal{L}_{B}+\mathcal{L}_{O},$

where

$\displaystyle\mathcal{L}_{A}:\
\dot{\mathbf{Y}}=\mathbf{U}\mathbf{U}^{T}\dot{\mathbf{A}},$ (15)
$\displaystyle\mathcal{L}_{B}:\
\dot{\mathbf{Y}}=-\mathbf{U}\mathbf{U}^{T}\dot{\mathbf{A}}\mathbf{V}\mathbf{V}^{T},$
$\displaystyle\mathcal{L}_{O}:\
\dot{\mathbf{Y}}=\dot{\mathbf{A}}\mathbf{V}\mathbf{V}^{T}.$

######  Theorem 1.

For $\mathbf{Y}$ defined in Eqn. (3) and $\dot{\mathbf{Y}}$ defined in Eqn.
(4), assume that the condition described in Eqn.(5) is satisfied. Then the
analytical solution to the sub-generator $\mathcal{L}_{A}$ described in Eqn.
(15) is

$\displaystyle\dot{\overbrace{\mathbf{S}\mathbf{V}^{T}}}$
$\displaystyle=\mathbf{U}^{T}\dot{\mathbf{A}},$
$\displaystyle\dot{\mathbf{U}}$ $\displaystyle=\mathbf{0}.$

The analytical solution to the sub-generator $\mathcal{L}_{B}$ is

$\displaystyle\dot{\mathbf{S}}=-\mathbf{U}^{T}\dot{\mathbf{A}}\mathbf{V},$
$\displaystyle\dot{\mathbf{U}}=\mathbf{0},$
$\displaystyle\dot{\mathbf{V}}=\mathbf{0}.$

And the analytical solution to the sub-generator $\mathcal{L}_{O}$ is

$\displaystyle\dot{\overbrace{\mathbf{U}\mathbf{S}}}$
$\displaystyle=\dot{\mathbf{A}}\mathbf{V},$ $\displaystyle\dot{\mathbf{V}}$
$\displaystyle=\mathbf{0}.$

Based on Theorem 1, we devise a novel method for Problem (2). The splitting
integration scheme also allows us to alternatively update $\mathbf{U}$,
$\mathbf{V}$ and $\mathbf{S}$, rather than directly update $\mathbf{Y}$. Owing
to the Markovian property of the Kolmogorov operator, different orders of sub-
generators $\mathcal{L}_{A}$, $\mathcal{L}_{B}$ and $\mathcal{L}_{O}$ are
equivalent Leimkuhler and Matthews (2013). In our work, we restrict our
interest on ‘OBA’ scheme: $\mathcal{L}_{O}+\mathcal{L}_{B}+\mathcal{L}_{A}$.

We call our method the splitting numerical integration method. The detail is
given in Algorithm 1. Here
$\sigma_{\text{min}}(\mathbf{V}_{i-1}^{T}\mathbf{V}_{i})$ measures the
distance between column spaces of $\mathbf{V}_{i-1}$ and $\mathbf{V}_{i}$, and
it is employed as the stopping criteria.

Compared with dynamical low rank approximation in Eqn. (9) which directly
updates $\mathbf{Y}$ and the stepsize must be small enough (approaching to 0),
our splitting numerical integration method updates the three components
($\mathbf{U}$, $\mathbf{S}$ and $\mathbf{V}$) in a “finer-grained” manner. For
instance, more specifically, in Step 7 of Algorithm 1, we adopt the fresh
$\mathbf{U}_{i}$ rather than $\mathbf{U}_{i-1}$. Notice that in Algorithm 1
the stepsize is implicitly set to 1. In summary, our approach is able to
address the issues mentioned in the end of Section 2.1.

###  3.3 Convergence Analysis

In this section we study convergence properties of our splitting numerical
integration algorithm.

######  Theorem 2.

Let $\\{\mathbf{Y}_{i}\colon i=0,1,\ldots\\}$ denote the sequence generated by
Algorithm 1. Then

$\displaystyle\|\mathbf{Y}_{i-1}-\mathbf{M}\|_{F}>\|\mathbf{Y}_{i}-\mathbf{M}\|_{F}.$

Furthermore, it can be also shown that our splitting numerical integration
method converges in terms of the subspace estimation.

######  Theorem 3.

Let $\\{\mathbf{U}_{i}\colon i=0,1,\ldots\\}$ be the sequence generated by
splitting numerical integration in Algorithm 1. Then the subspace error
$\|\mathbf{U}_{i}\mathbf{U}_{i}^{T}\mathbf{M}-\mathbf{M}\|_{F}$ decreases
monotonically until convergence. In particular, we have

$\|\mathbf{U}_{i-1}\mathbf{U}_{i-1}^{T}\mathbf{M}-\mathbf{M}\|_{F}>\|\mathbf{U}_{i}\mathbf{U}_{i}^{T}\mathbf{M}-\mathbf{M}\|_{F}.$

The similar results hold for $\\{\mathbf{V}_{i}\colon i=0,1,\ldots\\}$. That
is,

$\|\mathbf{M}\mathbf{V}_{i-1}\mathbf{V}_{i-1}^{T}-\mathbf{M}\|_{F}>\|\mathbf{M}\mathbf{V}_{i}\mathbf{V}_{i}^{T}-\mathbf{M}\|_{F}.$

Algorithm 1 splitting numerical integration

0: matrix $\mathbf{M}\in\mathbb{R}^{m{\times}n}$, target rank $r$, initial
value
$\mathbf{Y}_{0}=\mathbf{U}_{0}\mathbf{S}_{0}\mathbf{V}_{0}^{T}\in\mathcal{M}_{r}$,
tolerance for stopping criteria $\tau<1$, maximal iteration $T$.

0: truncated SVD of $\mathbf{M}$.

1: for $i=1:T$ do

2: Compute the gradient $\dot{\mathbf{A}}=\nabla f(X_{i-1})$. # O($mn$) flops

3: $\mathbf{Q}=\dot{\mathbf{A}}\mathbf{V}_{i-1}\ \ \ \in\mathbb{R}^{m\times
r}$ # O($mnr$) flops

4: $\mathbf{K}=\mathbf{U}_{i-1}\mathbf{S}_{i-1}+\mathbf{Q}\ \ \
\in\mathbb{R}^{m\times r}$ # O($mr^{2}$) flops

5: Perform QR-factorization to $\mathbf{K}$:
$[\mathbf{U}_{i},\mathbf{S}_{i-1+1/3}]=\text{QR}(\mathbf{K})$ # O($mr^{2}$)
flops

6: $\mathbf{S}_{i-1+2/3}=\mathbf{S}_{i-1+1/3}-\mathbf{U}^{T}_{i}\mathbf{Q}\ \
\ \in\mathbb{R}^{r\times r}$ # O($mr^{2}$) flops

7:
$\mathbf{L}=\mathbf{V}_{i-1}\mathbf{S}_{i-1+2/3}^{T}+\dot{\mathbf{A}}^{T}\mathbf{U}_{i}\
\ \ \in\mathbb{R}^{n\times r}$ # O($mnr$) flops

8: Perform QR-factorization to $\mathbf{L}$:
$[\mathbf{V}_{i},\mathbf{S}_{i}^{T}]=\text{QR}(\mathbf{L})$ # O($nr^{2}$)
flops

9: if ($\sigma_{\text{min}}(\mathbf{V}_{i-1}^{T}\mathbf{V}_{i})>\tau$) then

10: break.

11: end if # O($nr^{2}$) flops

12: end for

13: Perform SVD on the matrix
$\mathbf{S}_{T}=\mathbf{U}_{s}\mathbf{D}\mathbf{V}_{s}^{T}$, then we have that
$\mathbf{M}_{r}=(\mathbf{U}_{T}\mathbf{U}_{s})\mathbf{D}(\mathbf{V}_{T}\mathbf{V}_{s})^{T}$.

###  3.4 Application in low-rank matrix completion

The matrix completion problem is to recover a low-rank matrix from a few
observations of this matrix. In fixed-rank formulation of matrix completion,
we modify Problem (2) into

$\displaystyle\underset{\mathbf{Y}}{\arg\min}\frac{1}{2}\|P_{\Omega}({\mathbf{M}})-P_{\Omega}({\mathbf{Y}})\|_{F}^{2}$
$\displaystyle\text{s.t.}\ \mathbf{Y}\in\mathcal{M}_{r},$

where $\Omega$ represents the index of observations. Naturally, our splitting
numerical integration method applies to this scenario with only a modification
that the objective function becomes

$f_{1}(\mathbf{Y})=\frac{1}{2}\|P_{\Omega}({\mathbf{M}})-P_{\Omega}({\mathbf{Y}})\|_{F}^{2}.$
(16)

Then the convergence property can be extended into the partial observation
case. In the following theorem, we prove that the objective function is
monotonically decreasing until reaching the convergence condition.

######  Theorem 4.

Let $\\{\mathbf{Y}_{i}\colon i=0,1,\ldots\\}$ be the sequence generated by
splitting numerical integration in Algorithm 1 under any observation index
$\Omega$. Then splitting numerical integration decreases monotonically in the
objective function $f_{1}$ defined in Eqn. (16); that is,

$\displaystyle f_{1}(\mathbf{Y}_{i-1})>f_{1}(\mathbf{Y}_{i}).$

##  4 Empirical Evaluation

In this section, we conduct the empirical analysis of the splitting numerical
integration method. First we analyze the performance of splitting numerical
integration for low rank matrix approximation on simulated datasets. Then we
validate the performance of splitting numerical integration for low rank
matrix completion on a set of real data Hoang et al. (2019); Xiao et al.
(2020).

###  4.1 Low rank matrix approximation

We evaluate splitting numerical integration on low rank matrix approximation
with comparison with some popular baseline methods. The baseline methods
contain the power method Gu (2015); Fu and others (2021), randomized SVD
(RSVD) Halko et al. (2011) and dynamical low rank approximation Koch and
Lubich (2007). The primary goal is to illustrate the approximate accuracy on
three simulated datasets.

In particular, the three target matrices with different ranks are randomly
generated. The error is measured by
$\|\mathbf{U}\mathbf{S}\mathbf{V}-\mathbf{M}_{r}\|_{F}/\|\mathbf{M}_{r}\|_{F}$.
Since the complexity of these approaches are all $O(mnr)$, the runtime is
similar for all methods and not reported. Each trial is conducted three
independent times and average error are reported in Table 1. We observe that
our splitting numerical integration method owns obvious advantage over other
baseline method in accuracy. In addition, we use a special initialization by
letting the column spaces of $\mathbf{U}_{0}$ and $\mathbf{V}_{0}$ lie in the
orthogonal complements of $\mathbf{U}_{r}$ and $\mathbf{V}_{r}$, respectively.
The similar accuracy can be obtained, which means that our splitting numerical
integration method is insensitive to initialization.

Table 1: Performance of all methods on low rank matrix approximation, DLRA denotes dynamical low rank approximation and RSVD denotes randomized SVD. Note that “(20K,20K,20K)” corresponds to number of row, column and rank, respectively.  Method | A:(20K,20K,20K) | B:(20K,20K,2K) | C:(20K,20K,200)
---|---|---|---
DLRA | 2.8e-03 | 1.0e-03 | 1.2e-03
RSVD | 9.3e-02 | 6.0e-02 | 7.2e-02
power | 2.5e-08 | 4.7e-07 | 7e-08
SNI | 2.3e-11 | 2.1e-10 | 1.4e-10

###  4.2 Low rank matrix completion

Table 2: Results of recommendation systems measured in terms of the RMSE. ‘-’ represents the absence of results, which means that corresponding algorithm fails on this task due to memory or running time issue. Data set | Soft-Impute | ALS | GECO | LMaFit | RP | ScGrass | LRGoemCG | SNI
---|---|---|---|---|---|---|---|---
Movielens 100K | 0.9026 | 0.9696 | 0.9528 | 1.0821 | 0.9508 | 0.9502 | 0.9643 | 0.9501
Movielens 1M | 0.9127 | 0.9159 | 0.8601 | 0.8972 | 0.8590 | 0.8723 | 0.8934 | 0.8612
Movielens 10M | 0.8915 | 0.8726 | 0.8241 | 0.8921 | 0.8290 | 0.8991 | 0.8779 | 0.823
Netflix | 0.9356 | 0.9501 | 0.8738 | 0.9247 | 0.8601 | 0.9232 | 0.8723 | 0.8612
Yahoo Music | 24.77 | 24.59 | - | 26.43 | 23.93 | - | 24.09 | 22.86
Table 3: Running time (in seconds) of all methods on recommendation systems.  Data set | Soft-Impute | ALS | GECO | LMaFit | RP | ScGrass | LRGoemCG | SNI
---|---|---|---|---|---|---|---|---
Movielens 100K | 2.56 | 0.49 | 2.90 | 0.230 | 0.21 | 0.92 | 0.99 | 0.094
Movielens 1M | 22.81 | 5.61 | 176.11 | 1.412 | 1.00 | 50.11 | 15.23 | 0.94
Movielens 10M | 675.11 | 88.40 | $>10^{3}$ | 159.80 | 147.34 | $>10^{3}$ | 313.30 | 47.93
Netflix | $>5\times 10^{3}$ | 1189.47 | $>10^{4}$ | 345.00 | 744.43 | $>5\times 10^{3}$ | 3823.13 | 350.38
Yahoo Music | $>5\times 10^{4}$ | 8522.23 | - | 1239.56 | 1858.43 | - | 4043.22 | 236.32

We now conduct the empirical analysis of our splitting numerical integration
method for the low rank matrix completion (LRMC) problem. To show the
efficiency and effectiveness of our splitting numerical integration-LRMC, we
compare it with a bunch of baseline methods, including Soft-Impute Mazumder et
al. (2010), ALS (Soft-Impute Alternating Least Squares) Hastie et al. (2014),
GECO (Greedy Efficient Component Optimization) Shalev-Shwartz et al. (2011),
LMaFit (Low Rank Matrix Fitting) Wen et al. (2010), RP (Riemann Pursuit for
matrix recovery) Tan et al. (2014), ScGrass (Scaled Gradient on Grassmann
Manifold) Ngo and Saad (2012), and LRGeomCG (Low rank Geometric Conjugate
Gradient) Vandereycken (2013). The codes of all the methods can be available
online, e.g., Soft-Impute and ALS 111http://web.stanford.edu/ hastie/pub.htm,
GECO222http://www.cs.huji.ac.il/ shais/code/index.html,
LMaFit333http://lmafit.blogs.rice.edu/, RP and
LRGoemCG444http://www.tanmingkui.com/rp.html, and ScGrass555http://www-
users.cs.umn.edu/ thango/. These algorithms have been proved to be state-of-
the-art algorithms in low rank matrix completion.

We compare these methods on several popular recommendation systems. It is
worth mentioning that large-scale recommendation systems (say, Yahoo Music and
Netflix) are used to evaluate the scalability of our method. We use five
publicly available datasets: Movielens 100K, 1M, 10M, NetFlix, Yahoo Music
Track 1 to evaluate both the effectiveness and efficiency of our method.

Testing error in terms of RMSE (Root-Mean-Square Error) and computational
efficiency measured by running time are shown in Table 2 and Table 3,
respectively. From Table 2, we can observe that our method can achieve better
performance than most of the baseline methods in terms of RMSE Fu and others
(2020); Gao et al. (2019). Whilst Table 3 shows that our method can achieve
great speedup compared with almost all baseline methods under the same
setting. It is worth mentioning that in large scale tasks such as Yahoo Music
dataset, some results are not listed, which means that the corresponding
algorithm can not handle these cases in limited time or simply fail in these
situations.

## References

  * Drineas et al. [2006] Petros Drineas, Michael W Mahoney, and S Muthukrishnan.  Subspace sampling and relative-error matrix approximation: Column-based methods.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 316–326. Springer, 2006.
  * Eckart and Young [1936] Carl Eckart and Gale Young.  The approximation of one matrix by another of lower rank.  Psychometrika, 1(3):211–218, 1936.
  * Fu and others [2020] Tianfan Fu et al.  $\alpha$-mop: Molecule optimization with $\alpha$-divergence.  In 2020 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pages 240–244. IEEE, 2020.
  * Fu and others [2021] Tianfan Fu et al.  Probabilistic and dynamic molecule-disease interaction modeling for drug discovery.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pages 404–414, 2021.
  * Gao et al. [2019] Tian Gao, Tengfei Ma, and Jimeng Sun.  Pearl: Prototype learning via rule learning.  In Proceedings of the 10th ACM International Conference on Bioinformatics, Computational Biology and Health Informatics, pages 223–232, 2019.
  * Golub and Van Loan [2012] Gene H Golub and Charles F Van Loan.  Matrix computations, volume 3.  JHU Press, 2012.
  * Gu [2015] Ming Gu.  Subspace iteration randomization and singular value problems.  SIAM Journal on Scientific Computing, 37(3):A1139–A1173, 2015.
  * Hairer et al. [2006] Ernst Hairer, Marlis Hochbruck, Arieh Iserles, and Christian Lubich.  Geometric numerical integration.  Oberwolfach Reports, 3(1):805–882, 2006.
  * Halko et al. [2011] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM review, 53(2):217–288, 2011.
  * Hastie et al. [2014] Trevor Hastie, Rahul Mazumder, Jason Lee, and Reza Zadeh.  Matrix completion and low-rank svd via fast alternating least squares.  arXiv preprint arXiv:1410.2596, 2014.
  * Hoang et al. [2019] Trong Nghia Hoang, Cao Xiao, and Jimeng Sun.  Ddl: Deep dictionary learning for predictive phenotyping.  In IJCAI: proceedings of the conference, volume 2019, page 5857, 2019.
  * Huang and others [2020] Kexin Huang et al.  Deeppurpose: A deep learning library for drug-target interaction prediction.  Bioinformatics, 2020.
  * Jolliffe [2002] Ian Jolliffe.  Principal component analysis.  Wiley Online Library, 2002.
  * Koch and Lubich [2007] Othmar Koch and Christian Lubich.  Dynamical low-rank approximation.  SIAM Journal on Matrix Analysis and Applications, 29(2):434–454, 2007.
  * Leimkuhler and Matthews [2013] Benedict Leimkuhler and Charles Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Applied Mathematics Research eXpress, 2013(1):34–56, 2013.
  * Mazumder et al. [2010] Rahul Mazumder, Trevor Hastie, and Robert Tibshirani.  Spectral regularization algorithms for learning large incomplete matrices.  The Journal of Machine Learning Research, 11:2287–2322, 2010.
  * Muller et al. [2004] Neil Muller, Lourenço Magaia, and Ben M Herbst.  Singular value decomposition, eigenfaces, and 3d reconstructions.  SIAM review, 46(3):518–545, 2004.
  * Ngo and Saad [2012] Thanh Ngo and Yousef Saad.  Scaled gradients on grassmann manifolds for matrix completion.  In Advances in Neural Information Processing Systems, pages 1412–1420, 2012.
  * Shalev-Shwartz et al. [2011] Shai Shalev-Shwartz, Alon Gonen, and Ohad Shamir.  Large-scale convex minimization with a low-rank constraint.  arXiv preprint arXiv:1106.1622, 2011.
  * Tan et al. [2014] Mingkui Tan, Ivor W Tsang, Li Wang, Bart Vandereycken, and Sinno J Pan.  Riemannian pursuit for big matrix recovery.  In Proceedings of the 31st International Conference on Machine Learning (ICML-14), pages 1539–1547, 2014.
  * Vandereycken [2013] Bart Vandereycken.  Low-rank matrix completion by riemannian optimization.  SIAM Journal on Optimization, 23(2):1214–1236, 2013.
  * Wen et al. [2010] Z Wen, W Yin, and Y Zhang.  Low-rank factorization model for matrix completion by a non-linear successive over-relaxation algorithm.  Technical report, Rice CAAM Tech Report 10-07, University of Rice, 2010\.
  * Xiao et al. [2020] Cao Xiao, Jimeng Sun, and Xinhao Li.  MIMOSA: Multi-constraint molecule sampling for molecule optimization.  the Association for the Advancement of Artificial Intelligence, 2020\.