# Adaptive and Implicit Regularization for Matrix Completion

Zhemin Li
Department of Mathematics
National University of Defense Technology
Changsha 410008,Peoples Republic of China
<EMAIL_ADDRESS>
&Tao Sun
College of Computer
National University of Defense Technology
Changsha 410008,Peoples Republic of China
<EMAIL_ADDRESS>
&Hongxia Wang
Department of Mathematics
National University of Defense Technology
Changsha 410008,Peoples Republic of China
<EMAIL_ADDRESS>corresponding author
&Bao Wang
Department of Mathematics, Scientific Computing and Imaging Institute
University of Utah, Salt Lake City, UT, 84112, USA
<EMAIL_ADDRESS>

###### Abstract

The explicit low-rank regularization, e.g., nuclear norm regularization, has
been widely used in imaging sciences. However, it has been found that implicit
regularization outperforms explicit ones in various image processing tasks.
Another issue is that the fixed explicit regularization limits the
applicability to broad images since different images favor different features
captured by different explicit regularizations. As such, this paper proposes a
new adaptive and implicit low-rank regularization that captures the low-rank
prior dynamically from the training data. The core of our new adaptive and
implicit low-rank regularization is parameterizing the Laplacian matrix in the
Dirichlet energy-based regularization, which we call the regularization AIR.
Theoretically, we show that the adaptive regularization of AIR enhances the
implicit regularization and vanishes at the end of training. We validate AIR’s
effectiveness on various benchmark tasks, indicating that the AIR is
particularly favorable for the scenarios when the missing entries are non-
uniform. The code can be found at https://github.com/lizhemin15/AIR-Net.

_K_ eywords low-rank regularization, adaptive and implicit regularization,
matrix completion, deep learning

##  1 Introduction

Natural images usually lie in a high-dimensional ambient space, but with a
much lower intrinsic dimension [1], which has motivated many statistical
priors for image processing, including sparsity [2] and low-rank [3]. Other
priors, such as Total Variation (TV) and Dirichlet Energy (DE), also play
vital roles in image processing. There is no standard way to choose a proper
prior to a particular task, and people always depend on empirical experiences.

This paper is devoted to establishing a new adaptive regularization. The core
of our proposed adaptive regularization is parameterizing the Laplacian matrix
in the conventional DE regularization, and this regularization changes on the
fly as the underlying matrix gets updated. In particular, the proposed
adaptive regularizer estimates the prior matrix based on the historical
information during learning iterations for matrix completion. Compared to the
existing regularization schemes for matrix completion, the proposed
regularizer is adaptable to matrix completion problems from different
applications with different missing patterns. Furthermore, our adaptive
regularizer learns hyperparameters from data rather than manually tuning,
saving tremendous computational costs. Theoretically, we prove that solving
our proposed model with gradient descent can enhance the low-rank property of
the recovered matrices, c.f. Theorem 1, and the gradient descent iterations
will result in a non-trivial minimum, c.f. Theorem 2. Numerically, we verify
our proposed adaptive regularizer’s superiority over existing ones on various
benchmark matrix completion tasks.

We organize this paper as follows: we first introduce the preliminary for
readers unfamiliar with the regularization methods in Section 2. We will
integrate our proposed adaptive regularization scheme with DMF in Section 3.
Then we illustrate and theoretically prove two implicit regularization
properties of our proposed adaptive regularization scheme in Section 4. We
empirically validate the efficacy of our proposed adaptive regularization
scheme in Section 5, followed by concluding remarks. The missing details of
technical proofs are provided in the appendix.

##  2 Preliminary

###  2.1 Related works

####  2.1.1 Low-rank matrix completion

Modeling and learning the low-dimensional and low-rank structures of natural
images is an important and exciting research area in signal processing and
machine learning communities [4, 5, 6, 7, 8]. The low-rank structure is
fascinating and widely used prior to image processing. One of the major
obstacles in enforcing low-rank prior for image processing is that the exact
low-rank minimization is a discrete optimization problem and is NP-hard [9].
The computational challenge of direct low-rank minimization has motivated
several surrogate models of the low-rank regularization. In particular, the
authors of [10] relax the low-rank minimization to the nuclear norm
minimization, which can be further relaxed to convex optimization with linear
constraints. Nevertheless, the convex relaxation solves the low-rank
minimization problem while sacrificing the model’s accuracy significantly.
Indeed, the convex relaxation can introduce unacceptable errors, especially
for the problems that are insufficiently low-rank[11]. In response,
researchers have shifted their attention to the nonconvex models for low-rank
guarantees. For instance, the matrix factorization model has been proposed to
compensate for the model of nuclear norm relaxation [12, 13, 14, 15]. The
matrix factorization model can be formulated as follows:

Given a matrix ${\bm{X}}^{*}\in{\mathbb{R}}^{m\times n}$, we seek the
decomposition ${\bm{X}}={\bm{W}}^{[0]}\cdot{\bm{W}}^{[1]}$ and ensure that
${\bm{X}}\approx{\bm{X}}^{*}$, where ${\bm{W}}^{[0]}\in{\mathbb{R}}^{m\times
r}$ and ${\bm{W}}^{[1]}\in{\mathbb{R}}^{r\times n}$ with $r\leq\min\\{m,n\\}$.
The decomposition ${\bm{W}}^{[0]}\cdot{\bm{W}}^{[1]}$ ensures the rank of
${\bm{X}}$ to be capped by $r$, and such a decomposition can be numerically
computed by alternating the update of ${\bm{W}}^{[0]}$ and ${\bm{W}}^{[1]}$.
It is straightforward to generalize the above matrix decomposition to the
product of multiple matrices (a.k.a. multi-layer matrix factorization), i.e.,
${\bm{X}}=\prod_{\ell=0}^{L-1}{\bm{W}}^{[\ell]}$111For the sake of simplicity,
we denote the decomposition as
${\bm{X}}=\prod_{\ell=0}^{L-1}{\bm{W}}^{[\ell]}$., which is called the Deep
Matrix Factorization (DMF) [16]. Interestingly, the multi-layer matrix
factorization enjoys implicit low-rank regularization without any special
requirement on the initialization and constraint on $r$, which overcomes the
sensitivity on initialization in the two-layer matrix factorization [17].

####  2.1.2 Implicit regularization: deep matrix factorization

Due to the advantage of DMF over many existing matrix completion algorithms,
various efforts have been contributed to this area. The seminal research about
implicit regularization originates from [18], which considers the recovery of
a positive semi-definite matrix from symmetric measurements. Based on the
shallow matrix factorization model ${\bm{X}}={\bm{W}}^{[0]}{\bm{W}}^{[1]}$,
the authors of [18] prove that when ${\bm{W}}^{[0]},{\bm{W}}^{[1]}$ are
initialized to $\alpha{\bm{I}}$ with $\alpha>0$, their methods can find the
minimal nuclear norm when $\alpha\rightarrow 0$. In [16], Arora et al. prove a
similar result for arbitrary depth and asymmetric cases; however, they ignore
the depth effects. As the main contribution of [16], Arora et al. relax the
assumptions on the initialization and get rid of the assumption of symmetric
measurements. Therefore, we can use more general and convenient
initializations that factors in the factorization is balanced at the
initialization, i.e.,
$\left({\bm{W}}^{[j+1]}(0)\right)^{\top}{\bm{W}}^{[j+1]}(0)={\bm{W}}^{[j]}(0)\left({\bm{W}}^{[j]}(0)\right)^{\top}$
for $j=1,\ldots,L-1$. It is worth mentioning that the initialization adopted
in [18] also satisfy the balanced initialization condition. The theoretical
results in [16] are established by using the gradient dynamics of DMF. In
particular, the singular values of
${\bm{X}}(t)=\prod_{\ell=0}^{L-1}{\bm{W}}^{[\ell]}(t)$ evolve with different
speeds and thus induce the implicit low-rank regularization. More
interestingly, Arora et al. find that the implicit regularization becomes more
significant as $L$ increases [16].

Nevertheless, the theoretically-principled initialization methods mentioned
above are not popular in practice, instead, Gaussian initialization is
preferred for practical usage [19]. The bottleneck of the Gaussian
initialization is that it may not yield a low-rank bias starting from the
Neural Tangent Kernel (NTK) regime 222This phenomenon happens when DMF is
sufficiently wide, and the variance of the initialization is large enough.
[20, 21, 22, 23, 24, 25]. However, outside the NTK regime, there exists an
active regime where the dynamics of DMF are nonlinear and favor the low-rank
solution [26, 27, 28, 29, 30, 31, 32]. Jacot et al. further discuss the
effects of initialization for implicit low-rank regularization in [33]. As the
norm of initialized parameters goes to zero, Jacot et al. conjecture that the
trajectory of the gradient flow goes from one saddle point to another saddle
point, corresponding to matrices with increased ranks. Therefore, we only need
to initialize the parameters of DMF to be sufficiently small for the two-layer
matrix factorization model. Another advantage of DMF is that we do not need to
estimate $r$ for matrix factorization, instead we can simply choose
$r=\min\left\\{m,n\right\\}$. Empirically, DMF outperforms the nuclear norm
regularization and the two-layer matrix factorization models for low-rank
matrix completion. As shown in Figure 1 (a), DMF achieves remarkable
performance for exact low-rank matrix completion problems, and the performance
becomes better as the depth of DMF model increases. Moreover, DMF has been
applied to many other practical applications, including natural image
processing and recommendation systems [34, 35, 11].

|  |
---|---|---
(a) Random missing | (b) Random missing | (c) Random missing
|  |
(d)Patch missing | (e) Textural missing | (f) Random missing
Figure 1: Trajectories of optimizing DMF with different depths for different
matrix completion problems, where we set the mean squared error (MSE) on
observed elements to be less than $10^{-3}$ as the stop criterion. Obs_MSE and
Un_MSE stand for the MSE on observed and unobserved elements, respectively.
(a) DMF for completing a rank-5 random matrix of size $100\times 100$ with
$80\%$ random missing entries; (b) DMF for completing the gray scale Barbara
image, which has the size $240\times 240$ but $80\%$ of its randomly sampled
entries are missing; (c) DMF + TV for completing the gray scale Barbara image,
which has the size $240\times 240$ but $80\%$ of its randomly sampled entries
are missing; (d) DMF + TV for completing the gray scale Barbara image with
missing patch (see Figure 9 (c)), which has the size $240\times 240$; (e) DMF
+ TV for completing the gray scale Barbara image with missing textural (see
Figure 9 (b)), which has the size $240\times 240$; (f) DMF + TV for completing
the IC ([11]), which has the size $240\times 240$ but $80\%$ of its randomly
sampled entries are missing.

####  2.1.3 Total variation and Dirichlet energy regularization

Compared to solving exact low-rank matrix completion problems, directly
applying DMF for solving low-rank-related practical problems, including image
inpainting, often gives opposite results. As shown in Figure 1 (b), the deep
matrix factorization model performs inferior to the shallow models. By
investigating the optimization trajectory during the gradient descent
procedure, we notice that the major issue that causes deeper models to perform
worse than the shallow ones is because the underlying image inpainting problem
cannot be precisely formulated as a low-rank matrix completion problem.
Indeed, some detailed structures of the matrix correspond to small singular
values of the matrix, which the low-rank prior cannot entirely model. To
resolve the above issue, we need to enforce other priors to compensate for
subtleties in solving image inpainting problems and beyond.

A simple but valuable prior for natural image processing is piece-wise
smoothness, which can be induced using the Total Variation (TV) regularization
[36, 37, 38]. TV has been integrated with DMF, achieving promising recovery
results for matrix completion problems, especially when the rate of missing
entries is high [34]. As illustrated in Figure 1 (c), deep DMF models
outperform shallow DMF models when TV regularization is enforced. Though TV is
a natural prior for natural images, it is suboptimal or inappropriate for many
other practical applications, including recommendation systems through matrix
completion. Besides the TV regularization, the Dirichlet Energy (DE)-based
regularization, defined through the Laplacian matrix [11], has also been
employed for low-rank matrix completion. The DE regularization can describe
the similarity between rows and columns of a given matrix. Moreover, we can
view the DE regularization as a nonlocal extension of the TV regularization.

###  2.2 Do we need a new regularization: Fixed or adaptive

Given a matrix completion problem, we can integrate several existing explicit
regularizers with DMF or other matrix completion models to obtain some
successful results. However, the design or selection of explicit regularizers
is often task- and missing pattern-dependent. As shown in Figure 1 (d), (e),
and (f), although TV regularization performs quite well for random missing
patterns and natural images, the performance of using TV regularizers
dramatically degrades when they are applied for other missing patterns and
data. As such, we need to develop a new regularizer adaptable to different
missing patterns and different types of data.

###  2.3 Notations

We denote scalars by lower or upper case letters; vectors and matrices by
lower and upper case boldface letters, respectively. For a vector
${\bm{x}}=(x_{1},\cdots,x_{d})^{\top}\in\mathbb{R}^{d}$, we use
$\|{\bm{x}}\|:={(\sum_{i=1}^{d}|x_{i}|^{2})^{1/2}}$ to denote its
$\ell_{2}$-norm. We denote the matrix whose entries are all 1s as
$\mathbf{1}_{m,n}\in\mathbb{R}^{m\times n}$. For a matrix ${\bm{A}}$, we use
${\bm{A}}^{\top}$, ${\bm{A}}^{-1}$, and $\|{\bm{A}}\|$ to denote its
transpose, inverse, and spectral norm, respectively. Given a matrix
${\bm{A}}$, we denote ${\bm{A}}_{ij}$ as its $(i,j)$-th entry, and denote its
$k$-th column and $k$-th row as ${\bm{A}}_{:,k}$ and ${\bm{A}}_{k,:}$,
respectively. For a function
$f({\bm{x}}):\mathbb{R}^{d}\rightarrow\mathbb{R}$, we denote $\nabla
f({\bm{x}})$ as its gradient.

##  3 Adaptive Laplacian Regularization

Besides low-rank, self-similarity is another widely used prior for solving
inverse problems, and the self-similarity is usually measured by the DE. A
particular type of DE is the TV [36], which describes the piece-wise self-
similarity within a given image. To adapt the DE regularization to different
missing patterns and different types of data, we raise its degree of freedom
by introducing a learnable DE regularization. In particular, we propose the
following model

$\textrm{minimize}_{{\bm{X}},{\bm{W}}_{i}}\left\\{{\mathcal{L}}_{\textrm{all}}:=\Bigg{(}{\mathcal{L}}_{{\mathbb{Y}}}\left({\bm{Y}},{\mathcal{A}}({\bm{X}})\right)+\sum_{i=1}^{N}\lambda_{i}\cdot{\mathcal{R}}_{{\bm{W}}_{i}}\left({\mathcal{T}}_{i}\left({\bm{X}}\right)\right)\Bigg{)}\right\\},$
(1)

where the linear operator ${\mathcal{A}}(\cdot):\mathbb{R}^{m\times
n}\mapsto\mathbb{R}^{o}$, ${\bm{Y}}\in{\mathbb{Y}}\subseteq\mathbb{R}^{o}$ is
a measure vector and ${\mathcal{L}}_{\mathbb{Y}}$ is a distance metric on
${\mathbb{Y}}$. Different from other regularization models for matrix
completion, here _${\bm{X}}\in\mathbb{R}^{m\times n}$ is parameterized by
products of matrix which tends to be low-rank implicitly_, and
${\mathcal{R}}_{{\bm{W}}_{i}}$ is an adaptive regularizer and $\lambda_{i}\geq
0$ is a constant. ${\mathcal{T}}_{i}(\cdot)$ is a linear or nonlinear
transformation. See Section 3.2 for details. A particular case of (1) for
matrix completion is given in Section 3.3 below.

###  3.1 Self-similarity and Dirichlet energy

Given a matrix ${\bm{X}}\in{\mathbb{R}}^{m\times n}$, DE associated with the
weighted adjacency matrix ${\bm{A}}\in{\mathbb{R}}^{m\times m}$, along rows of
${\bm{X}}$, is a classical approach to encode the self-similarity prior of the
matrix ${\bm{X}}$. The adjacency matrix ${\bm{A}}$ measures the similarity
between rows of ${\bm{X}}$, and ${\bm{A}}_{ij}$ is bigger if rows $i$ and $j$
of ${\bm{X}}$ are more similar. The corresponding Laplacian matrix of
${\bm{A}}$ is ${\bm{L}}:={\bm{D}}-{\bm{A}}$, where ${\bm{D}}$ is the degree
matrix with ${\bm{D}}_{ii}=\sum_{j=1}^{n}{\bm{A}}_{ij}$ and ${\bm{D}}_{ij}=0$
if $i\neq j$. Then we can mathematically formulate the DE as follows

$\text{tr}({\bm{X}}^{\top}{\bm{L}}{\bm{X}})=\sum_{1\leq i,j\leq
m}{\bm{A}}_{ij}\left\|{\bm{X}}_{i,:}-{\bm{X}}_{j,:}\right\|^{2}.$

It is evident that when DE is minimized, rows ${\bm{X}}_{i,:}$ and
${\bm{X}}_{j,:}$ become closer with bigger ${\bm{A}}_{ij}$.

There are two major obstacles to using DE in applications: (i) ${\bm{L}}$ is
unknown for an incomplete matrix, and (ii) DE only encodes the similarity
between rows of ${\bm{X}}$; other similarities such as block similarity cannot
be captured. To resolve the above two issues, we parameterize ${\bm{L}}$ as a
function of ${\bm{X}}$ and learn it during completing ${\bm{X}}$.

###  3.2 Learned Dirichlet energy as an adaptive regularizer

Directly parameterizing ${\bm{L}}$ without enforcing any particular structure
and then optimizing it by minimizing
$\text{tr}({\bm{X}}^{\top}{\bm{L}}{\bm{X}})$ might result in all entries of
${\bm{L}}$ being very negative. To remedy this problem, we first recall a few
essential properties of ${\bm{L}}$ as follows:

  1. 1.

The Laplacian matrix ${\bm{L}}$ can be generated from the corresponding
adjacency matrix ${\bm{A}}$, that is ${\bm{L}}={\bm{D}}-{\bm{A}}$, where
${\bm{D}}$ is the degree matrix.

  2. 2.

${\bm{A}}$ is a symmetric matrix and ${\bm{A}}_{ij}\geq 0$ for all $i,j$.

One of the most straightforward approaches is to model ${\bm{A}}$ as the
covariance matrix, e.g., ${\bm{A}}:={\bm{X}}^{\top}{\bm{X}}$, of the rows of
${\bm{X}}$. However, this model may not capture sufficient similarities among
rows of ${\bm{A}}$. To this end, the self-attention mechanism [39] has been
proposed to model similarities among different rows. However, the learning
self-attention mechanism usually requires a large amount of training data,
which is not the case for the matrix completion problem. Therefore, we cannot
leverage the off-the-shelf self-attention mechanism for matrix completion.

We propose a new parameterization of the adjacency matrix that is similar to
but simpler than the self-attention mechanism. Our new parametrization can
capture row similarities of a given matrix beyond the covariance. In
particular, we parameterize ${\bm{A}}$ as follows

$\left\\{\begin{array}[]{l}{\bm{A}}=\frac{\exp({\bm{W}}+{\bm{W}}^{\top})}{\mathbf{1}_{m}^{\top}\exp({\bm{W}})\mathbf{1}_{m}}\\\
{\bm{L}}=\left({\bm{A}}\cdot\mathbf{1}_{m\times
m}\right)\odot\mathbf{I}_{m}-{\bm{A}}\end{array}\right.,$ (2)

where $\mathbf{1}_{m}$ is an $m$-dimensional vector whose entries are all 1s,
$\mathbf{I}_{m}$ is the $m\times m$ identity matrix, $\mathbf{1}_{m\times m}$
is the $m\times m$ matrix whose entries are all 1s and $\odot$ is Hadamard
product. ${\bm{W}}\in{\mathbb{R}}^{m\times m}$ has the same size as ${\bm{A}}$
and $\exp(\cdot)$ denotes the element-wise exponential. Note that the above
parameterization of ${\bm{A}}$ guarantees ${\bm{A}}$ to be symmetric, and all
its entries are positive. To search for the optimal ${\bm{A}}$, we need to
update ${\bm{W}}$ start from a certain initialization. It is clear that the
minimum value of DE is 0, i.e.,
$\text{tr}({\bm{X}}^{\top}{\bm{L}}{\bm{X}})=0$, which is obtained when
${\bm{A}}={\bm{L}}={\bf 0}$ and this is known as the trivial solution. Our
parameterization (2) directly avoids the above trivial solution since
${\bm{A}}_{ij}>0$. In Theorem 2, we will show that ${\bm{L}}$ will converge to
a non-trivial minimum under the gradient descent dynamics using the
parameterization in (2).

To capture more types of self-similarity, we further replace ${\bm{X}}$ with
${\mathcal{T}}_{i}\left({\bm{X}}\right)$, where ${\mathcal{T}}_{i}(\cdot)$ is
a linear or nonlinear transformation that will be specified in particular
applications. And we can generalize the above adaptive regularizer as follows

${\mathcal{R}}_{{\bm{W}}_{i}}\left({\mathcal{T}}_{i}\left({\bm{X}}\right)\right)=\mathrm{tr}\left(\left[{\mathcal{T}}_{i}\left({\bm{X}}\right)\right]^{\top}{\bm{L}}_{i}\left({\bm{W}}_{i}\right){\mathcal{T}}_{i}\left({\bm{X}}\right)\right),i=1,2,\ldots,N,$
(3)

where ${\bm{L}}_{i}\in{\mathbb{R}}^{m_{i}\times m_{i}}$ is parameterized by
${\bm{W}}_{i}\in{\mathbb{R}}^{m_{i}\times m_{i}}$ following (2).

${\mathcal{T}}_{i}:{\mathbb{R}}^{m\times n}\mapsto{\mathbb{R}}^{m_{i}\times
n_{i}}$ transforms ${\bm{X}}$ into another domain, making the adaptive
regularization able to capture different similarities in data. The common
choice of ${\mathcal{T}}_{i}(\cdot)$ can be
${\mathcal{T}}_{i}\left({\bm{X}}\right)={\bm{X}}$ and
${\mathcal{T}}_{i}\left({\bm{X}}\right)={\bm{X}}^{\top}$, which capture the
row and column correlation, respectively. Another particularly interesting
case is
${\mathcal{T}}_{i}({\bm{X}})=\left[\textbf{vec}\left(\text{block}({\bm{X}})\right)_{1},\ldots,\textbf{vec}\left(\text{block}({\bm{X}})\right)_{m_{i}}\right]^{\top}$,
where
$\textbf{vec}\left(\text{block}({\bm{X}})\right)_{j}\in{\mathbb{R}}^{n_{i}}$
is the vectorization of $j$-th block of ${\bm{X}}$ (we divide ${\bm{X}}$ into
$m_{i}$ blocks), then the similarity among blocks can be modeled. Apart from
these handcraft transformations to model specific structured correlations of
the image, we can also parameterize ${\mathcal{T}}_{i}$ with a neural network
to learn more broad classes of correlations from data. However, we noticed
that learning ${\mathcal{T}}_{i}$ using a neural network is a very challenging
task, which we leave as future work. This work focuses on handcrafted
transformations ${\mathcal{T}}_{i}$ that capture similarities among rows and
columns of the matrix ${\bm{X}}$. Another interesting result is that
${\mathcal{R}}_{{\bm{W}}_{i}}$ vanishes at the end of the training, which
avoids over-fitting without early stopping. We show that
${\mathcal{R}}_{{\bm{W}}_{i}}$ vanishes at the end of training both
theoretically (see Theorem 2) and numerically (see Section 5).

###  3.3 Adaptive regularizer for matrix completion

This subsection considers applying our proposed adaptive regularizer for
matrix completion. In particular, we consider two transformations:
${\mathcal{T}}_{1}\left({\bm{X}}\right)={\bm{X}}$ and
${\mathcal{T}}_{2}\left({\bm{X}}\right)={\bm{X}}^{\top}$ to capture row and
column similarities of ${\bm{X}}$, respectively. We denote the corresponding
subscript of row and column as $r$ and $c$, respectively. And the DMF model
with the regularization in (1) becomes

$\displaystyle\textrm{minimize}{\bm{{\bm{W}}^{[l]},{\bm{W}}_{r},{\bm{W}}_{c}}}$
$\displaystyle\left\\{{\mathcal{L}}:=\Bigg{(}\frac{1}{2}\left\|{\bm{Y}}-{\mathcal{A}}\left(\prod_{l=0}^{L-1}{\bm{W}}^{[l]}\right)\right\|_{F}^{2}\right.$
(4)
$\displaystyle\left.+\lambda_{r}\cdot{\mathcal{R}}_{{\bm{W}}_{r}}\left(\prod_{l=0}^{L-1}{\bm{W}}^{[l]}\right)+\lambda_{c}\cdot{\mathcal{R}}_{{\bm{W}}_{c}}\left(\left(\prod_{l=0}^{L-1}{\bm{W}}^{[l]}\right)^{\top}\right)\Bigg{)}\right\\},$

for $l=0,1,\cdots,L-1.$ We name (4) as Adaptive and Implicit Regularization
(AIR). The parameters in AIR are updated by using gradient descent or its
variants. We stop the iteration when
$\left|{\mathcal{R}}_{{\bm{W}}_{r}(T+1)}-{\mathcal{R}}_{{\bm{W}}_{r}(T)}\right|<\delta$
and
$\left|{\mathcal{R}}_{{\bm{W}}_{c}(T+1)}-{\mathcal{R}}_{{\bm{W}}_{c}(T)}\right|<\delta$.
The recovered matrix is
$\hat{{\bm{X}}}(T)={\bm{W}}^{[L-1]}(T)\cdots{\bm{W}}^{[0]}(T)$, which tends to
be low-rank implicitly [16]. With the implicit low-rank property, we do not
need to estimate the shared dimension of ${\bm{W}}^{[l]}$ in advance (more
details can be found in Appendix C). The framework of algorithm is shown in
Figure 2.

Some works that combine implicit and explicit regularization can be regarded
as special cases of (1). For instance, both TV [36] and DE [40] can be
considered as fixed ${\bm{L}}$ in (1). Therefore, our proposed framework in
(1) contains DMF+TV [34] and DMF+DE [11].

Figure 2: We give an example of AIR when completing a $3\times 5$ matrix. The
${\mathcal{A}}$ and ${\bm{Y}}$ are inputs of our algorithm. In this case, we
factorized ${\bm{X}}$ into three matrices. The ${\bm{W}}_{r}$ and
${\bm{W}}_{c}$ are parameters of AIR. In each iteration, we calculate the loss
function based on these parameters and update them with an optimization
algorithm such as gradient descent.

##  4 Theoretical Analysis of AIR

In this section, we will analyze the theoretical properties of solving problem
(4) using gradient descent. In particular, we will show that (a) AIR enhances
the implicit low-rank of DMF (see Theorem 1 below); and (b) the adaptive
regularization will converge to a minimum under gradient descent, and the
resulting model captures the intrinsic structure of data flexibly. Although we
focus on matrix completion problems, the following theoretical analyses apply
to general inverse problems. For the ease of notation, as ${\mathcal{A}}$ and
${\bm{Y}}$ are fixed during optimization, we denote
${\mathcal{L}}_{{\mathbb{Y}}}\left({\bm{Y}},{\mathcal{A}}({\bm{X}})\right)$ as
${\mathcal{L}}_{{\mathbb{Y}}}\left({\bm{X}}\right)$ below.

###  4.1 Adaptive regularizer has implicit low-rank regularization

We show that the proposed adaptive regularizer introduces implicit low-rank
regularization. Below we denote the Laplacian matrices parameterized by
${\bm{W}}_{r}$ and ${\bm{W}}_{c}$, following (2), as ${\bm{L}}_{r}$ and
${\bm{L}}_{c}$, respectively.

######  Theorem 1.

Consider the following dynamics with initial data satisfying the balanced
initialization in Assumption 1 (see Appendix A for details)

$\displaystyle\dot{{\bm{W}}}^{[l]}(t)=-\frac{\partial}{\partial{\bm{W}}^{[l]}}{\mathcal{L}}({\bm{X}}(t)),\quad
t\geq 0,\quad l=0,\ldots,L-1,$

where
${\mathcal{L}}\left({\bm{X}}\right):={\mathcal{L}}_{{\mathbb{Y}}}({\bm{X}})+\lambda_{r}\cdot{\mathcal{R}}_{{\bm{W}}_{r}}\left({\bm{X}}\right)+\lambda_{c}\cdot{\mathcal{R}}_{{\bm{W}}_{c}}\left({\bm{X}}\right)$.
Then for $k=1,2,\ldots$, we have

$\displaystyle\dot{\sigma}_{k}(t)=$
$\displaystyle-L\left(\sigma_{k}^{2}(t)\right)^{1-\frac{1}{L}}\left\langle\nabla_{{\bm{W}}}{\mathcal{L}}_{{\mathbb{Y}}}({\bm{X}}(t)),{\bm{U}}_{k,:}(t)\left[{\bm{V}}_{k,:}(t)\right]^{\top}\right\rangle-2L\left(\sigma_{k}^{2}(t)\right)^{\frac{3}{2}-\frac{1}{L}}\gamma_{k}(t),$
(5)

where ${\bm{U}}(t){\bm{S}}(t){\bm{V}}(t)^{\top}$ is the SVD of the matrix
${\bm{X}}\left(t\right)$, $\sigma_{k}(t)$ is the $k$-th singular value in
descending order of ${\bm{X}}(t)$,
${\bm{X}}=\prod_{l=0}^{L-1}{\bm{W}}^{[l]}=\sum\limits_{s}\sigma_{s}{\bm{U}}_{s,:}{\bm{V}}_{s,:}^{\top},\gamma_{k}(t)={\bm{U}}_{k,:}^{\top}{\bm{L}}_{r}{\bm{U}}_{k,:}+{\bm{V}}_{k,:}^{\top}{\bm{L}}_{c}{\bm{V}}_{k,:}\geq
0$.

###### Proof.

By direct calculation, we have

$\displaystyle\nabla_{{\bm{X}}}\left(\lambda_{r}{\mathcal{R}}_{{\bm{W}}_{r}}+\lambda_{c}{\mathcal{R}}_{{\bm{W}}_{c}}\right)$
$\displaystyle=\frac{\partial\mathrm{tr}\left(\lambda_{r}{\bm{X}}^{\top}{\bm{L}}_{r}{\bm{X}}+\lambda_{c}{\bm{X}}{\bm{L}}_{c}{\bm{X}}^{\top}\right)}{\partial{\bm{X}}}=2\lambda_{r}{\bm{L}}_{r}{\bm{X}}+2\lambda_{c}{\bm{X}}{\bm{L}}_{c}$
$\displaystyle=2\lambda_{r}{\bm{L}}_{r}\sum\limits_{s}\sigma_{s}{\bm{U}}_{s,:}{\bm{V}}_{s,:}^{\top}+2\lambda_{c}\sum\limits_{s}\sigma_{s}{\bm{U}}_{s,:}{\bm{V}}_{s,:}^{\top}{\bm{L}}_{c}.$

Note that

$\langle{\bm{V}}_{s,:},{\bm{V}}_{s^{\prime},:}\rangle=\langle{\bm{U}}_{s,:},{\bm{U}}_{s^{\prime},:}\rangle=\delta_{ss^{\prime}}=\left\\{\begin{array}[]{cc}1,&s=s^{\prime}\\\
0,&s\neq s^{\prime}\\\ \end{array}\right.,$

Therefore

$\displaystyle{\bm{U}}_{k,:}^{\top}\left(\nabla_{{\bm{X}}}\left(\lambda_{r}{\mathcal{R}}_{{\bm{W}}_{r}}+\lambda_{c}{\mathcal{R}}_{{\bm{W}}_{c}}\right)\right){\bm{V}}_{k,:}$
$\displaystyle=2\sigma_{k}(\lambda_{r}{\bm{U}}_{k,:}^{\top}{\bm{L}}_{r}{\bm{U}}_{k,:}+\lambda_{c}{\bm{V}}_{k,:}^{\top}{\bm{L}}_{c}{\bm{V}}_{k,:})=2\sigma_{k}\gamma_{k}(t),$

where the term
$\gamma_{k}(t)=2\sigma_{k}(\lambda_{r}{\bm{U}}_{k,:}^{\top}{\bm{L}}_{r}{\bm{U}}_{k,:}+\lambda_{c}{\bm{V}}_{k,:}^{\top}{\bm{L}}_{c}{\bm{V}}_{k,:})\geq
0$. Furthermore, according to (7) in the appendix, we have

${\bm{U}}_{k,:}^{\top}\nabla_{{\bm{X}}}{\mathcal{L}}{\bm{V}}_{k,:}=-L\left(\sigma_{k}^{2}(t)\right)^{1-\frac{1}{L}}\left\langle\nabla_{{\bm{X}}}{\mathcal{L}}({\bm{X}}(t)),{\bm{U}}_{k,:}(t)\left[{\bm{V}}_{k,:}(t)\right]^{\top}\right\rangle.$

Combined the above results with
${\mathcal{L}}({\bm{X}})={\mathcal{L}}_{{\mathbb{Y}}}({\bm{X}})+\lambda_{r}{\mathcal{R}}_{{\bm{W}}_{r}}+\lambda_{c}{\mathcal{R}}_{{\bm{W}}_{c}}$,
we have

$\dot{\sigma}_{k}(t)=-L\left(\sigma_{k}^{2}(t)\right)^{1-\frac{1}{L}}\left\langle\nabla_{{\bm{X}}}{\mathcal{L}}_{{\mathbb{Y}}}({\bm{X}}(t)),{\bm{U}}_{k,:}(t)\left[{\bm{V}}_{k,:}(t)\right]^{\top}\right\rangle-2L\left(\sigma_{k}^{2}(t)\right)^{\frac{3}{2}-\frac{1}{L}}\gamma_{k}(t),$

which completes the proof. ∎

Compared with the result of the vanilla DMF whose order of $\sigma_{k}(t)$ is
$2-\frac{2}{L}$, Theorem 1 demonstrates that AIR’s $\sigma_{k}(t)$ has a
higher dynamics order $3-\frac{2}{L}$. Note that the adaptive regularizer
keeps $\gamma_{k}(t)\geq 0$. This way, a faster convergence rate gap appears
between different singular values $\sigma_{r}$ than the vanilla DMF.
Therefore, AIR enhances the implicit low-rank regularization over DMF.

###  4.2 The auto vanishing property of the adaptive regularizer

We suppose the matrix ${\bm{X}}$ is fixed and then study the convergence of
${\mathcal{R}}_{{\mathcal{T}}_{i}({\bm{W}}_{i})}$ based on the dynamics of
${\bm{L}}_{i}$. Theorem 2 shows that
${\mathcal{R}}_{{\mathcal{T}}_{i}({\bm{W}}_{i})}$ vanishes at the end and thus
prevents AIR from suffering over-fitting. In the process of proof of the
following, we can replace the ${\mathcal{T}}_{i},{\bm{W}}_{i}$ with
${\mathcal{T}},{\bm{W}}$ in the following proof.

######  Theorem 2.

Consider the following gradient flow model,

$\dot{{\bm{W}}}(t)=-\nabla_{{\bm{W}}(t)}{\mathcal{R}}_{{\bm{W}}(t)}\left({\mathcal{T}}\left({\bm{X}}(t)\right)\right)=-\nabla_{{\bm{W}}(t)}\mathrm{tr}\left({\mathcal{T}}\left({\bm{X}}(t)\right)^{\top}{\bm{L}}\left({\bm{W}}(t)\right){\mathcal{T}}\left({\bm{X}}(t)\right)\right),$

where $\|{\mathcal{T}}\left({\bm{X}}\right)_{k,:}\|_{F}^{2}=1$ and
${\mathcal{T}}\left({\bm{X}}\right)_{kl}>0$. If we initialize
${\bm{W}}\left(0\right)=\varepsilon\mathbf{1}_{m_{i}\times
m_{i}},\forall\varepsilon\in\mathbb{R}$, then ${\bm{W}}\left(t\right)$ will
keep symmetric during the optimization procedure. Furthermore, we have the
following element-wise convergence

$\left|{\bm{L}}_{kl}(t)-{\bm{L}}_{kl}^{*}\right|\leq\left\\{\begin{array}[]{cc}2\exp(-Dt)/\gamma,&(k,l)\in{\mathbb{S}}_{1}\\\
\exp\left(-Dt\right),&(k,l)\in{\mathbb{S}}_{2}\\\
2\left(m-1\right)\exp(-Dt)/\gamma,&k=l\end{array}\right.,$

where

${\mathbb{S}}_{1}=\left\\{(k,l)\mid k\neq
l,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}\neq{\mathcal{T}}\left({\bm{X}}\right)_{l,:}\right\\},\quad{\mathbb{S}}_{2}=\left\\{(k,l)\mid
k\neq
l,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}={\mathcal{T}}\left({\bm{X}}\right)_{l,:}\right\\},$

and

${\bm{L}}_{kl}^{*}=\left\\{\begin{array}[]{cc}0,&(k,l)\in{\mathbb{S}}_{1}\\\
\gamma,&(k,l)\in{\mathbb{S}}_{2}\\\ -\sum_{l^{\prime}=1,l^{\prime}\neq
l}^{m}{\bm{L}}_{kl^{\prime}}^{*},&k=l\end{array}\right..$

${\bm{L}}_{kl}(t)$ is the $(k,l)$-th element of ${\bm{L}}(t)$,
$\mathbf{1}_{m\times m}$ is a matrix whose entries are all 1s. $\gamma<1,D\geq
0$ are constant defined in Appendix B and $D$ equals zero if and only if
${\mathcal{T}}\left({\bm{X}}\right)=\textbf{1}_{m_{i}\times n_{i}}$.

###### Proof.

See Appendix B. ∎

The assumptions in Theorem 2 can be satisfied by applying appropriate
normalization and linear transformation on ${\bm{X}}$. Theorem 2 gives the
limit and convergence rate of ${\bm{L}}(t)$. In particular, the limit
satisfies ${\bm{L}}^{*}_{kl}=0$ unless
${\mathcal{T}}\left({\bm{X}}\right)_{k,:}={\mathcal{T}}\left({\bm{X}}\right)_{l,:}$
or $k=l$. That is, $\mathbf{L}^{*}_{i,j}$ reflects if the $i$-th and $j$-th
rows of the matrix $\mathcal{T}\left(\mathbf{X}\right)$ are the same or not.
$\mathbf{L}^{*}_{i,j}\neq 0$ if and only if the $i$-th and $j$-th rows of
$\mathcal{T}\left(\mathbf{X}\right)$ are identical. ${\bm{L}}_{kl}(t)$
converges faster for those $(k,l)\in{\mathbb{S}}_{2}$ than that in
${\mathbb{S}}_{1}$. In other words, the adaptive regularizer first captures
structural similarities in the matrix ${\bm{X}}$. The different convergence
rates for different $(k,l)$ indicate AIR generates a multi-scale similarity,
which will be discussed in Section 5.1. Moreover,
${\mathcal{R}}_{{\bm{W}}}(t)$ converges to 0 as shown in the following
Corollary. The adaptive regularization will vanish at the end of training and
not cause over-fitting.

######  Corollary 1.

In the setting of Theorem 2, we further have

$0\leq{\mathcal{R}}_{{\bm{W}}}(t)\leq
2\left(m-1\right)m\exp(-Dt)/\gamma,i=1,2,\ldots,N.$

###### Proof.

Direct computation gives

$\displaystyle{\mathcal{R}}_{{\bm{W}}}(t)$
$\displaystyle=\sum_{k,l}{\bm{L}}_{kl}(t)\left\|{\mathcal{T}}({\bm{X}})_{k,:}-{\mathcal{T}}({\bm{X}})_{l,:}\right\|_{2}^{2}$
$\displaystyle\leq 2\sum_{k\neq l}{\bm{L}}_{kl}(t)\leq
2m\left(m-1\right)\exp(-Dt)/\gamma,$

which concludes the proof. ∎

We have shown that AIR can enhance implicit low-rank regularization and avoid
over-fitting. In the next section, we will validate these theoretical results
and the effectiveness of AIR numerically.

##  5 Experimental Results

In this section, we numerically validate the following adaptive properties of
AIR: (a) Laplacian matrices ${\bm{L}}_{r}$ and ${\bm{L}}_{c}$ capture the
structural similarity in data from large scale to small scale; (b) the
capability of capturing structural similarities at all scales is crucial to
the success of matrix completion, which confirms the importance of our
proposed adaptive regularizers; (c) AIR is adaptive to different data, avoids
overfitting, and achieves remarkable performance; (d) AIR behaves like the
momentum that has been widely used in optimization [41, 42, 43, 44, 45] and
neural network architecture design [46, 47].

Data types and sampling patterns. We consider completion of three types of
matrices: gray-scale image, user-movie rating matrix [11], and drug-target
interaction (DTI) data [11]. In particular, we consider three benchmark images
of size $240\times 240$ [48], including Baboon, Barbara, and Cameraman. We
consider the Syn-Netflix dataset for the user-movie rating matrix, which has a
size of $150\times 200$. The DTI data describes the interaction of Ion
channels (IC) and G protein-coupled receptor (GPCR), which have the shape of
$210\times 204$ and $223\times 95$, respectively [11, 49]. Moreover, we study
matrix completion with three different missing patterns: random missing, patch
missing, and textural missing, see Figure 9 for an illustration. The random
missing rate varies in different experiments, and the default missing rate is
$30\%$.

Parameter settings. We set
$\lambda_{r}=\lambda_{c}={({\bm{Y}}_{\max}-{\bm{Y}}_{\min})}/{(mn)}$ to ensure
both fidelity and regularization have similar order of magnitude, where
${\bm{Y}}_{\max}$ and ${\bm{Y}}_{\min}$ are the maximum and the minimum entry
of ${\bm{Y}}$, respectively. $\delta$ is a threshold, which has the default
value of ${mn}/{1000}$. All the parameters in AIR are initialized with
Gaussian initialization of zero mean and variance $10^{-5}$. We use Adam [50]
to train the AIR.

###  5.1 Multi-scale similarity captured by adaptive regularizer

In this subsection, we will validate Theorems 1 and 2. We focus on using AIR
for inpainting the corrupted Baboon image and completing the Syn-Netflix
matrix. Figure 3 and Figure 4 plot the heatmaps of Laplacian matrices
${\bm{L}}_{r}(t)$ and ${\bm{L}}_{c}(t)$ for Baboon and Syn-Netflix
experiments, respectively. These heatmaps show what AIR learns during
training.

Figure 3: (a): first and second rows depict the Baboon image and its rotation. (b)-(d): first/second row shows the heatmap of ${\bm{L}}_{r}$/${\bm{L}}_{c}$ at different $t$s. A darker color indicates a stronger similarity captured by the adaptive regularizer. The $(i,j)$-th element in the heatmap of ${\bm{L}}_{r}(t)$ has a darker color than the $(i,j^{\prime})$-th element indicates that the $i$-th row is more related to $j$-th row compared with $j^{\prime}$-th row.  |  |  |
---|---|---|---
(a) Users relationship | (b) ${\bm{L}}_{r}(t=5000)$ | (c) ${\bm{L}}_{r}(t=7000)$ | (d) ${\bm{L}}_{r}(t=10000)$
|  |  |
(e) Items relationship | (f) ${\bm{L}}_{c}(t=5000)$ | (g) ${\bm{L}}_{c}(t=7000)$ | (h) ${\bm{L}}_{c}(t=10000)$
Figure 4: The two panels in the first column shows the exact correlation among
rows and columns, respectively. The other three columns are the Laplacian
matrix learned by AIR at different iterations.

As shown in Figure 3, both Laplacian matrices ${\bm{L}}_{r}(t)$ and
${\bm{L}}_{c}(t)$ appear many blocks during the initial training stage, e.g.,
$t=4000$. We consider a few boxed blocks of the matrix ${\bm{L}}_{c}$ in
Figure 3 (b, bottom); the values in these blocks reflect the strong similarity
of the corresponding highlighted patches of the original Baboon image, which
echos our intuition. Moreover, the slight difference between those selected
columns will be future captured by the adaptive regularizer as the training
goes, e.g., at $t=7000$. As the training goes further, the regularization
gradually vanishes. We see that the regularization becomes rather weak at
$t=10000$. Syn-Netflix experiments also support the above conclusions,
${\bm{L}}_{r}$ and ${\bm{L}}_{c}$ tend to the Users/Items relationship (in
Figure 4 (a, e)). Different from the case in Figure 3, the ${\bm{L}}_{r}(t)$
and ${\bm{L}}_{c}(t)$ in Figure 4 retain the value when the rating matrix in
Figure 8 (a) corresponding rows or columns are identical.

These results confirm that AIR captures the similarity from large to small.
This raises an interesting question: _does there exist an optimal $t^{*}$ such
that ${\bm{L}}_{r}(t^{*})$ and ${\bm{L}}_{c}(t^{*})$ are precisely captured by
AIR? If yes, can we use these fixed optimal ${\bm{L}}_{r}(t^{*})$ and
${\bm{L}}_{c}(t^{*})$ for AIR and further train ${\bm{X}}$ to obtain even
better matrix completion?_

Next, we experimentally show that Laplacian matrices ${\bm{L}}_{r}$ and
${\bm{L}}_{c}$ learned by AIR are crucial for matrix completion. We fix
${\bm{L}}_{r}$ and ${\bm{L}}_{c}$ that were learned at a specific training
step for AIR. In particular, we compare AIR with adaptive ${\bm{L}}_{r}$ and
${\bm{L}}_{c}$, and the variants of AIR with both Laplacian matrices fixed as
that learned at $t=4000$, $7000$, and $9000$, respectively. We use the same
hyperparameters for training AIR with fixed Laplacian matrices as that used
for training AIR. To quantitatively compare the performance of AIR with
different settings, we consider the following Normalized Mean Absolute Error
(NMAE) to measure the gap between the recovered and exact elements at
unobserved locations,

$\mathrm{NMAE}=\frac{\left\|\bar{{\mathcal{A}}}({\bm{X}})-\bar{{\bm{Y}}}\right\|_{F}^{2}}{(mn-o)\left({\bm{Y}}^{*}_{\max}-{\bm{Y}}^{*}_{\min}\right)},$

where $\bar{{\bm{Y}}}$ is the unobserved elements, $\bar{{\mathcal{A}}}$ maps
${\bm{X}}$ to the corresponding unobserved location, and
${\bm{Y}}^{*}=\left[{\bm{Y}},\bar{{\bm{Y}}}\right]^{\top}$ is the ground
truth.

We contrast the vanilla AIR and AIR with fixed Laplacian matrices for Baboon
image inpainting. Figure 5 shows how the NMAE changes during training. AIR,
which updates the regularization during training, achieves the best
performance for all missing patterns. Fixing Laplacian matrices can accelerate
the convergence of training AIR at the beginning. While the optimal time
$t^{*}$ is unknown for AIR with fixed Laplacian matrices, the learned
Laplacian matrices work without estimating $t^{*}$.

|  |
---|---|---
(a) Patch | (b) Random | (c) Fixed
Figure 5: Contrasting the adaptive regularizer with the fixed regularizer for
Baboon image recovery. We consider inpainting the Baboon image with three
types of missing pixels: (a) patch missing, (b) randomly missing 30% pixels,
and (c) texture missing. The blue lines plot the NMAE during training the
vanilla AIR. The remaining three lines in each figure indicate replacing
${\bm{L}}_{r}$ and ${\bm{L}}_{c}$ with ${\bm{L}}_{r}(t)$ and ${\bm{L}}_{c}(t)$
at $t=4000$, $7000$, and $9000$, respectively.

###  5.2 AIR behaves like a momentum

This subsection makes a heuristic and direct connection between AIR and the
momentum method. AIR converges to the vanilla DMF model when the adaptive
regularization vanishes, and Corollary 1 shows that the adaptive
regularization vanishes exponentially fast. We need to distinguish AIR and DMF
according to the optimization trajectory. We compare three models (DMF,
DMF+TV, and DMF+AIR, i.e., AIR) in Figure 6, and we see that the three models
perform dramatically different near the convergence. At the beginning of
training, the observed and unobserved MSEs of the three models drop similarly.
When the observed MSE becomes smaller, the model learns details in observed
elements. The unobserved MSE increased during the observed MSE decrease in the
vanilla DMF and DMF+TV cases. Our proposed AIR keeps the decaying trend for
both observed and unobserved MSEs.

Looking back into the training of AIR, we see that the update of
${\bm{X}}(t+1)$ involves both ${\bm{X}}(t)$ and ${\bm{L}}(t)$, and the update
of ${\bm{L}}(t)$ depends on ${\bm{X}}(t-1)$. To understand the training
dynamics, we consider the following simplified model

$\textrm{minimize}_{{\bm{X}},{\bm{L}}}\left\\{{\mathcal{L}}={\mathcal{L}}_{\mathbb{Y}}+\text{tr}\left({\bm{X}}^{\top}{\bm{L}}{\bm{X}}\right)\right\\},$

and we have

$\nabla_{{\bm{X}}(t)}{\mathcal{L}}=\nabla_{{\bm{X}}(t)}{\mathcal{L}}_{\mathbb{Y}}+2{\bm{L}}(t){\bm{X}}(t),$

where ${\bm{L}}(t)$ is the function of $\left\\{{\bm{X}}(t_{0})\mid
t_{0}<t\right\\}$. Therefore, every iteration step of AIR leverages all the
learned previous information $\left\\{{\bm{X}}(t_{0})\mid t_{0}<t\right\\}$.
While the update of both vanilla DMF and DMF+TV only depends on ${\bm{X}}(t)$.
From this viewpoint, AIR shares a similar spirit as the momentum method, which
leverages history to improve performance.

|  |
---|---|---
(a) Patch | (b) Random | (c) Fixed
Figure 6: Trajectories of training DMF model without regularization, with TV,
and with AIR regularization for inpainting the Baboon image with three types
of missing pixels: (a) patch missing, (b) randomly missing 30% pixels, and (c)
texture missing. The dot point indicates $X(0)$ and the diamond shape of
different color indicate $X(10000)$ for different models.

###  5.3 Adaptive regularizer performance on varied data and missing pattern

We apply AIR for matrix completion on three data types with different missing
patterns. We will also compare AIR with a few baseline algorithms. Figure 7
plots the raw matrices of Syn-Netflix and DTI datasets.

Baseline algorithms. In the following experiments, we compare AIR with several
popular matrix completion algorithms, including KNN [51], SVD [52], PNMC [53],
DMF [16], and RDMF [34]. For Syn-Netflix-related experiments, we replace RDMF
with DMF+DE [11], which is a better fit for that task.

Avoid Over-fitting. Figure 8 shows the evolution of training NMAE of DMF and
AIR. AIR avoids over-fitting and performs better on all three data types with
various missing patterns compared with vanilla DMF.

|  |
---|---|---
(a) Syn-Netflix | (b) IC | (c) GPCR
Figure 7: Matrix completion for several benchmark datasets, including Syn-Netflix, IC, and GPCR. Deeper color in Syn-Netflix heatmap indicate a higher rate from a user to a given item. IC and GPCR are both drug-target interaction (DTI) [49] matrix, and the entries of the matrix represent the interaction between drug and target with deeper color represents stronger indication. We consider binary interaction for both IC and GPCR.  |  |
---|---|---
(a) Random Cameraman | (b) Textural Cameraman | (c) Patch Cameraman
|  |
(d) Random Syn-Netflix | (e) Random IC | (f) Random GPCR
Figure 8: Evolution of NMAE during training DMF and AIR (DMF with AIR regularization). The first row plots the results of Cameraman with (a) random missing, (b) textural missing, and (c) patch missing. The second row plots three other benchmarks with random missing: (d) Syn-Netflix, (e) IC, and (f) GPCR  Table 1:  NMAE values of different algorithms for solving different matrix completion problems with different missing patterns. NA indicates that the method is not suitable for that task. The hyper-parameters of models and algorithms keep consistent with the original paper. Data | Missing | KNN[51] | SVD[52] | PNMC[53] | DMF[16] | RDMF[34] | DMF+DE[11] | AIR
---|---|---|---|---|---|---|---|---
Barbara | $30\%$ | 0.083 | 0.0621 | 0.0622 | 0.0613 | 0.0494 | NA | 0.0471
Patch | 0.1563 | 0.2324 | 0.2055 | 0.7664 | 0.3025 | NA | 0.1195
Texture | 0.0712 | 0.1331 | 0.1100 | 0.3885 | 0.1864 | NA | 0.0692
Baboon | $30\%$ | 0.0831 | 0.1631 | 0.0965 | 0.2134 | 0.0926 | NA | 0.0814
Patch | 0.1195 | 0.1571 | 0.1722 | 0.8133 | 0.2111 | NA | 0.1316
Texture | 0.1237 | 0.1815 | 0.1488 | 0.5835 | 0.2818 | NA | 0.1208
Syn-Netflix | $70\%$ | 0.0032 | 0.0376 | NA | 0.0003 | NA | 0.0008 | 0.0002
$75\%$ | 0.0046 | 0.0378 | NA | 0.0004 | NA | 0.0009 | 0.0003
$80\%$ | 0.0092 | 0.0414 | NA | 0.0014 | NA | 0.0012 | 0.0007
IC | $20\%$ | 0.0169 | 0.0547 | NA | 0.0773 | NA | 0.0151 | 0.0134
GPCR | $20\%$ | 0.0409 | 0.0565 | NA | 0.1513 | NA | 0.0245 | 0.0271

Adaptive to data. Table 1 lists the NMAEs of matrix completion using AIR and
several benchmark algorithms for different data with different missing
patterns. We see that, in general, AIR performs the best and even better than
the well-calibrated algorithms for some particular datasets. We further plot
the recovered images in Figure 9; which shows that AIR consistently gives
appealing results visually. Some baseline algorithms also perform well for
some specific missing patterns. For instance, RDMF performs well for the
random missing case but poorly for the other missing patterns. PNMC can
complete the missing patch well but does not work well when the texture is
missing. Overall, AIR achieves decent results visually and quantitatively for
all data under all missing patterns.

Figure 9: Contrasting KNN [51], SVD [52], PNMC [53], DMF [16], RDMF [34], and
our proposed AIR on the Babara image inpainting with three types of missing
data respectively. The hyper-parameters of benchmark models and algorithms are
adopted from the original paper.

###  5.4 Contrasting computational time of AIR with DMF

The proposed AIR has several advantages, including avoiding over-fitting and
adapting to data. However, AIR has more parameters than vanilla DMF, requiring
more computational time in each iteration. In this subsection, we compare the
complexity of AIR with two benchmark algorithms, namely, DMF and RDMF.

We focus on square matrices of the size $m\times m$. The loss function of
$L$-layer DMF is given by
$\frac{1}{2}\left\|{\mathcal{A}}\left(\prod_{l=0}^{L-1}{\bm{W}}^{[l]}\right)-{\bm{Y}}^{*}\right\|_{F}^{2}$,
where ${\bm{W}}^{[l]}\in\mathbb{R}^{m\times m}$ for $l=0,\ldots,L-1$ and
${\bm{X}}=\prod_{l=0}^{L-1}{\bm{W}}^{[l]}$ is the recovered matrix. Therefore
the computational complexity of computing the loss function of DMF is
${\mathcal{O}}((L-1)m^{3})$. As for RDMF, its loss function involves the
calculation of an additional TV term whose computational complexity is
${\mathcal{O}}(m^{2})$. Therefore, RDMF is only slightly more expensive than
DMF. The AIR’s loss function calculates two additional adaptive regularization
terms $tr({\bm{X}}{\bm{L}}_{r}{\bm{X}}^{T})$ and
$tr({\bm{X}}^{T}{\bm{L}}_{c}{\bm{X}})$, with computational complexity being
${\mathcal{O}}(2m^{3})$.

The previous analysis shows that the computational complexity of each model is
dominated by the size of ${\bm{X}}$ ($m$) and the depth of the model ($L$). We
numerically verify the previous computational complexity analysis on the
Barbara image with different $m$ and $L$. In particular, we iterate each model
for $10000$ iterations and report the computational time in Table 2. We
conduct our experiments on Tesla V100 GPUs and report the GPU time consumption
of three kinds of DMF-based methods, where the reported GPU times are averaged
over ten independent runs. The numerical results show that:

  1. 1.

For any fixed $m$ and $L$, AIR takes more computational time than RDMF, and
DMF is computationally the cheapest among the three models.

  2. 2.

For any fixed $m$, as $L$ increases, the ratio between the computational time
of AIR and DMF decays. This is because the computational complexity of DMF
increases as DMF becomes deeper. Nevertheless, the computational complexity of
both TV and adaptive regularization is independent of the depth $L$.

  3. 3.

For any fixed depth $L$, as $m$ increases, the ratio between the computational
time of RDMF and DMF decays, while the ratio between the computational time of
AIR and DMF keeps near a constant. This result echoes our computational
complexity analysis. In particular, the computational complexity ratio between
RDMF and DMF is ${\mathcal{O}}(1+\frac{1}{(L-1)m})$ and the ratio between AIR
and DMF is ${\mathcal{O}}(1+\frac{2}{(L-1)})$.

Although AIR takes more computational time than the baseline algorithms, it
achieves significantly better accuracy than the baselines, which is valuable
in many applications.

Table 2: Contrasting computational time of DMF, RDMF, and AIR with different $m,L$. RDMF/DMF represents the ratio of the computational time of RDMF and DMF. AIR/DMF represents the ratio of the computational time of AIR and DMF. Unit: second. $m$ | $L$ | DMF[16] | RDMF[34] | RDMF/DMF | AIR | AIR/DMF
---|---|---|---|---|---|---
100 | 2 | 31.55 | 44.24 | 1.40 | 59.53 | 1.89
3 | 33.93 | 46.81 | 1.38 | 61.53 | 1.81
4 | 36.52 | 49.16 | 1.35 | 61.59 | 1.69
170 | 2 | 35.15 | 47.67 | 1.36 | 65.81 | 1.87
3 | 37.97 | 51.47 | 1.36 | 67.97 | 1.79
4 | 41.38 | 55.11 | 1.33 | 70.86 | 1.71
240 | 2 | 43.57 | 53.36 | 1.22 | 82.21 | 1.89
3 | 46.50 | 59.89 | 1.29 | 83.32 | 1.79
4 | 52.24 | 61.89 | 1.18 | 88.34 | 1.69

##  6 Conclusion

This paper proposes that AIR solve matrix completion problems without knowing
the prior in advance. AIR parameterizes the Laplacian matrix in DE and can
adaptively learn the regularization according to different data at different
training steps. In addition, we demonstrate that AIR can avoid the over-
fitting issue. AIR is a generic framework for solving inverse problems that
simultaneously encode the self-similarity and low-rank prior. Some interesting
properties of the proposed AIR deserve further theoretical analysis. Such as
the momentum phenomenon, which calls for a more detailed dynamic analysis.

##  Appendix A A Brief Review of DMF Algorithm and Theory

######  Assumption 1.

Factor matrices are balanced at the initialization, i.e.,

$\left[{{\bm{W}}^{[l+1]}}(0)\right]^{\top}{\bm{W}}^{[l+1]}(0)={\bm{W}}^{[l]}(0)\left[{{\bm{W}}^{[l]}}(0)\right]^{\top},\quad
l=0,\ldots,L-2.$

Under this assumption, Arora et al. have studied the gradient flow of the non-
regularized risk function ${\mathcal{L}}_{{\mathbb{Y}}}$, which is governed by
the following differential equation

$\dot{{\bm{W}}}^{[l]}(t)=-\frac{\partial}{\partial{\bm{W}}^{[l]}}{\mathcal{L}}_{{\mathbb{Y}}}\left({\bm{X}}(t)\right),\quad
t\geq 0,\quad l=0,\ldots,L-1,$ (6)

where the empirical risk ${\mathcal{L}}_{{\mathbb{Y}}}$ can be any analytic
function of ${\bm{X}}(t)$. According to the analyticity of
${\mathcal{L}}_{{\mathbb{Y}}}$, ${\bm{X}}(t)$ has the following singular value
decomposition

${\bm{X}}(t)={\bm{U}}(t){\bm{S}}(t)\left[{\bm{V}}(t)\right]^{\top},$

where
${\bm{U}}(t)\in{\mathbb{R}}^{m,\min\left\\{m,n\right\\}},{\bm{S}}(t)\in{\mathbb{R}}^{\min\left\\{m,n\right\\},\min\left\\{m,n\right\\}}$,
and ${\bm{V}}(t)\in{\mathbb{R}}^{\min\left\\{m,n\right\\},n}$ are analytic
functions of $t$; and for every $t$, the matrices ${\bm{U}}(t)$ and
${\bm{V}}(t)$ have orthonormal columns, while ${\bm{S}}(t)$ is diagonal. We
denote the diagonal entries of ${\bm{S}}(t)$ by $\sigma_{1}(t),\ldots,$
$\sigma_{\min\left\\{m,n\right\\}}(t)$, which are the signed singular values
of ${\bm{X}}(t)$. The columns of ${\bm{U}}(t)$ and ${\bm{V}}(t)$, denoted by
${\bm{U}}_{1}(t),\ldots,{\bm{U}}_{\min\left\\{m,n\right\\}}(t)$ and
${\bm{V}}_{1}(t),\ldots,{\bm{V}}_{\min\left\\{m,n\right\\}}(t),$ are the
corresponding left and right singular vectors, respectively. We have the
following result about the evolutionary dynamics of the singular values of
${\bm{X}}(t)$.

######  Proposition 1 ([16, Theorem 3]).

Consider the dynamics of (6) with initial data satisfying Assumption 1. Then
the dynamics of signed singular values $\sigma_{k}(t)$ of matrix ${\bm{X}}(t)$
is governed by the following equation

$\displaystyle\dot{\sigma}_{k}(t)=-L\left(\sigma_{k}^{2}(t)\right)^{1-\frac{1}{L}}\left\langle\nabla_{{\bm{X}}}{\mathcal{L}}_{{\mathbb{Y}}}({\bm{X}}(t)),{\bm{U}}_{k,:}(t)\left[{\bm{V}}_{k,:}(t)\right]^{\top}\right\rangle,\
k=1,\ldots,\min\left\\{m,n\right\\}.$ (7)

If the matrix factorization is non-degenerate, i.e., the factorization has
depth $L\geq 2,$ the singular values need not be signed (we may assume
$\sigma_{k}(t)\geq 0$ for all $t$).

Arora et al. claim that the term
$\left(\sigma_{k}^{2}(t)\right)^{1-\frac{1}{L}}$ enhances the dynamics of
large singular values while diminishing that of small ones. The
enhancement/diminishment becomes more significant as $L$ increases [16]. As
the ${\mathcal{L}}_{{\mathbb{Y}}}$ can be any analyticity empirical risk
function, we can replace ${\mathcal{L}}_{{\mathbb{Y}}}$ with arbitrary
analyticity function ${\mathcal{L}}$ in 7.

##  Appendix B Proof of Theorem 2

######  Proposition 2.

We have

$\nabla_{{\bm{W}}}\left({\mathcal{R}}_{{\bm{W}}}\left({\bm{X}}\right)\right)=2{\bm{C}}\odot{\bm{A}}-2\mathrm{tr}\left({\bm{C}}{\bm{A}}^{\prime}\right){\bm{A}}^{\prime},$

where
${\bm{A}}^{\prime}=\frac{\exp({\bm{W}}^{\top})}{\mathbf{1}_{m}^{\top}\exp({\bm{W}})\mathbf{1}_{m}}$,
${\bm{A}}={\bm{A}}^{\prime}+{{\bm{A}}^{\prime}}^{\top}$, and

${\bm{C}}=\left({\mathcal{T}}\left({\bm{X}}\right){{\mathcal{T}}\left({\bm{X}}\right)}^{\top}\odot{\bm{I}}_{m}\right)\mathbf{1}_{m\times
m}-{\mathcal{T}}\left({\bm{X}}\right){{\mathcal{T}}\left({\bm{X}}\right)}^{\top}.$

###### Proof.

We denote ${\bm{X}}={\mathcal{T}}\left({\bm{X}}\right)\in{\mathbb{R}}^{m\times
n}$, then we consider the differential of
$\mathrm{tr}\left({\bm{X}}^{\top}{\bm{L}}{\bm{X}}\right)$, i.e.,
$d\left[\mathrm{tr}\left({\bm{X}}^{\top}{\bm{L}}{\bm{X}}\right)\right]$, note
that

$\displaystyle
d\left[\mathrm{tr}\left({\bm{X}}^{\top}{\bm{L}}{\bm{X}}\right)\right]=\mathrm{tr}\left[d\left({\bm{L}}{\bm{X}}{\bm{X}}^{\top}\right)\right]$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left[\left(d{\bm{A}}\mathbf{1}_{m\times
m}\right)\odot
I_{m}{\bm{X}}{\bm{X}}^{\top}-d{\bm{A}}{\bm{X}}{\bm{X}}^{\top}\right]$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left[{\bm{X}}{\bm{X}}^{\top}\left({\bm{I}}_{m}\odot\left(d{\bm{A}}\mathbf{1}_{m\times
m}\right)\right)-{\bm{X}}{\bm{X}}^{\top}d{\bm{A}}\right]$ $\displaystyle=$
$\displaystyle\mathrm{tr}\left[\left({\bm{X}}{\bm{X}}^{\top}\odot{\bm{I}}_{m}\right)^{\top}d{\bm{A}}\mathbf{1}_{m\times
m}-{\bm{X}}{\bm{X}}^{\top}d{\bm{A}}\right]$ $\displaystyle=$
$\displaystyle\mathrm{tr}\left[\left(\left({\bm{X}}{\bm{X}}^{\top}\odot{\bm{I}}_{m}\right)\mathbf{1}_{m\times
m}\right)^{\top}d{\bm{A}}-{\bm{X}}{\bm{X}}^{\top}d{\bm{A}}\right]$
$\displaystyle=$ $\displaystyle\mathrm{tr}\left[\left(\mathbf{1}_{m\times
m}\left({\bm{X}}{\bm{X}}^{\top}\odot{\bm{I}}_{m}\right)-{\bm{X}}{\bm{X}}^{\top}\right)d{\bm{A}}\right]$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left[\left(\left({\bm{X}}{\bm{X}}^{\top}\odot{\bm{I}}_{m}\right)\mathbf{1}_{m\times
m}-{\bm{X}}{\bm{X}}^{\top}\right)d{\bm{A}}\right].$

We denote
${\bm{C}}=\left({\bm{X}}{\bm{X}}^{\top}\odot{\bm{I}}_{m}\right)\mathbf{1}_{m\times
m}-{\bm{X}}{\bm{X}}^{\top},S=\mathbf{1}_{m}^{\top}\exp\left({\bm{W}}\right)\mathbf{1}_{m}$,
then

$\displaystyle
d\left[\mathrm{tr}\left({\bm{X}}^{\top}{\bm{L}}{\bm{X}}\right)\right]=\mathrm{tr}\left({\bm{C}}d{\bm{A}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{S^{2}}\mathrm{tr}\left[{\bm{C}}\left(S\exp({\bm{W}}+{\bm{W}}^{\top})\odot
d({\bm{W}}+{\bm{W}}^{\top})\right)\right.$
$\displaystyle\left.-{\bm{C}}\left(\mathbf{1}_{m}^{\top}\left(\exp({\bm{W}})\odot
d{\bm{W}}\right)\mathbf{1}_{m}\right)\exp({\bm{W}}+{\bm{W}}^{\top})\right]$
$\displaystyle=$ $\displaystyle\mathrm{tr}\left[{\bm{C}}\left({\bm{A}}\odot
d\left({\bm{W}}+{\bm{W}}^{\top}\right)\right)\right]$
$\displaystyle-\frac{1}{S^{2}}\mathrm{tr}\left[\mathbf{1}_{m\times
m}\left(\exp({\bm{W}})\odot
d{\bm{W}}\right)\right]\mathrm{tr}\left[{\bm{C}}\exp\left({\bm{W}}+{\bm{W}}^{\top}\right)\right]$
$\displaystyle=$ $\displaystyle\mathrm{tr}\left[{\bm{C}}\left({\bm{A}}\odot
d\left({\bm{W}}+{\bm{W}}^{\top}\right)\right)\right]-\frac{1}{S}\mathrm{tr}\left(\exp({\bm{W}}^{T})d{\bm{W}}\right)\mathrm{tr}\left({\bm{C}}{\bm{A}}\right)$
$\displaystyle=$ $\displaystyle\mathrm{tr}\left[{\bm{C}}\left({\bm{A}}\odot
d\left({\bm{W}}+{\bm{W}}^{\top}\right)\right)\right]-\mathrm{tr}\left[\mathrm{tr}\left({\bm{C}}{\bm{A}}\right){\bm{A}}^{\prime}d{\bm{W}}\right]$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left[\left(\left({\bm{C}}^{\top}\odot{\bm{A}}\right)^{\top}+{\bm{C}}^{\top}\odot{\bm{A}}-\mathrm{tr}\left({\bm{C}}{\bm{A}}\right){\bm{A}}^{\prime}\right)d{\bm{W}}\right].$

Therefore,

$\displaystyle\nabla_{{\bm{W}}}\mathrm{tr}\left({\bm{X}}^{\top}{\bm{L}}_{i}{\bm{X}}\right)$
$\displaystyle=\left({\bm{C}}^{\top}\odot{\bm{A}}\right)^{\top}+{\bm{C}}^{\top}\odot{\bm{A}}-\mathrm{tr}\left({\bm{C}}{\bm{A}}\right){\bm{A}}^{\prime}$
$\displaystyle=2{\bm{C}}\odot{\bm{A}}-2\mathrm{tr}\left({\bm{C}}{\bm{A}}^{\prime}\right){\bm{A}}^{\prime}.$

Notice that
${\bm{X}}={\mathcal{T}}\left({\bm{X}}\right)\in{\mathbb{R}}^{m\times n}$, the
proposition is proved. ∎

###### Proof.

(Proof of Theorem 2) We rewrite the gradient in Proposition 2 into the
following element wise formulation:

$\dot{{\bm{W}}}_{kl}\left(t\right)=\left(2Ca(t)-4{\bm{C}}_{kl}\right){\bm{A}}_{kl}^{\prime}(t),$

where $Ca(t)=\mathrm{tr}\left({\bm{C}}(t){\bm{A}}^{\prime}(t)\right)$ and the
sub-index denotes the corresponding entry of the matrix.

With the assumption that
$\Big{\|}{\mathcal{T}}\left({\bm{X}}\right)_{k,:}\Big{\|}_{F}^{2}=1$ and
${\mathcal{T}}\left({\bm{X}}\right)_{kl}>0$, we have

$0\leq\left({\mathcal{T}}\left({\bm{X}}\right){\mathcal{T}}\left({\bm{X}}\right)^{\top}\right)_{kl}\leq
1.$

Therefore
${\bm{C}}_{kl}=\left({\mathcal{T}}\left({\bm{X}}\right){\mathcal{T}}\left({\bm{X}}\right)^{\top}\right)_{kk}-\left({\mathcal{T}}\left({\bm{X}}\right){\mathcal{T}}\left({\bm{X}}\right)^{\top}\right)_{kl}\geq
0$ and ${\bm{C}}_{kk}=0$ as
${\bm{A}}_{kl}^{\prime}=S^{-1}\exp({{\bm{W}}}_{kl})>0$, thus we have

$Ca(t)=\mathrm{tr}\left({\bm{C}}(t){\bm{A}}^{\prime}(t)\right)=\sum_{k=1,l=1}^{m}C_{kl}{\bm{A}}_{kl}^{\prime}(t)\geq
0.$

Denote ${\bm{C}}_{\hat{k}\hat{l}}\in\min\limits_{k,l}{\bm{C}}_{kl}$, we have
${\bm{C}}_{\hat{k}\hat{l}}\leq{\bm{C}}_{kl}$ and then consider

$\displaystyle\dot{{\bm{W}}}_{\hat{k}\hat{l}}(t)-\dot{{\bm{W}}}_{kl}(t)=2Ca(t)\left({\bm{A}}_{\hat{k}\hat{l}}^{\prime}(t)-{\bm{A}}_{kl}^{\prime}(t)\right)-4\left({\bm{C}}_{\hat{k}\hat{l}}{\bm{A}}_{\hat{k}\hat{l}}^{\prime}(t)-{\bm{C}}_{kl}{\bm{A}}_{kl}^{\prime}(t)\right).$

As we initialize ${\bm{W}}\left(0\right)=\varepsilon\mathbf{1}_{m\times m}$,
thus ${\bm{A}}_{kl}^{\prime}(0)=\frac{1}{m^{2}}$, for $\forall k,l$.
Therefore,

$\dot{{\bm{W}}}_{\hat{k}\hat{l}}(0)-\dot{{\bm{W}}}_{kl}(0)=-4\left({\bm{C}}_{\hat{k}\hat{l}}-{\bm{C}}_{kl}\right){\bm{A}}_{kl}^{\prime}(0)=-\frac{4}{m^{2}}\left({\bm{C}}_{\hat{k}\hat{l}}-{\bm{C}}_{kl}\right)\geq
0.$

Then we have ${\bm{W}}_{\hat{k}\hat{l}}(t)\geq{\bm{W}}_{kl}(t)$ and
${\bm{A}}_{\hat{k}\hat{l}}^{\prime}(t)\geq{\bm{A}}_{kl}^{\prime}(t)$, the
equality holds if and only if $t=0$ or
${\bm{C}}_{\hat{k}\hat{l}}={\bm{C}}_{kl}$. Furthermore,
$\dot{{\bm{W}}}_{\hat{k}\hat{l}}(t)-\dot{{\bm{W}}}_{kl}(t)\geq-4\left({\bm{C}}_{\hat{k}\hat{l}}-{\bm{C}}_{kl}\right){\bm{A}}_{kl}^{\prime}(0)$,
then
${\bm{W}}_{\hat{k}\hat{l}}(t)-{\bm{W}}_{kl}(t)\geq{\bm{D}}_{\hat{k},\hat{l},k,l}t$,
where
${\bm{D}}_{\hat{k},\hat{l},k,l}=-4\left({\bm{C}}_{\hat{k}\hat{l}}-{\bm{C}}_{kl}\right){\bm{A}}_{kl}^{\prime}(0)\geq
0$. Next, we consider

$\displaystyle{\bm{A}}_{\hat{l}\hat{k}}^{\prime}\left(t\right)$
$\displaystyle=\frac{\exp({\bm{W}}_{\hat{k}\hat{l}})}{\sum\limits_{k,l}\exp({\bm{W}}_{kl})}$
$\displaystyle=\frac{1}{\sum\limits_{k,l}\exp({\bm{W}}_{kl}-{\bm{W}}_{\hat{k}\hat{l}})}\geq\frac{1}{\sum\limits_{k,l}\exp\left(-{\bm{D}}_{\hat{k},\hat{l},k,l}t\right)}.$

As ${\bm{C}}_{kk}=0$ and ${\bm{C}}_{kl}\geq 0$, therefore
${\bm{C}}_{kk}\in\min\limits_{kl}{\bm{C}}_{kl}$=0. It is not difficult to show
that ${\bm{C}}_{\hat{k}\hat{l}}=0$ if and only if
${{\mathcal{T}}\left({\bm{X}}\right)}_{\hat{k},:}={{\mathcal{T}}\left({\bm{X}}\right)}_{\hat{l},:}$.
We further simplify our notations by introducing the following notations

${\mathbb{S}}_{1}=\left\\{(k,l)\mid k\neq
l,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}\neq{\mathcal{T}}\left({\bm{X}}\right)_{l,:}\right\\},{\mathbb{S}}_{2}=\left\\{(k,l)\mid
k\neq
l,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}={\mathcal{T}}\left({\bm{X}}\right)_{l,:}\right\\},$

${\mathbb{S}}_{3}=\left\\{(k,l)\mid k=l\right\\}.$

If we denote
${\mathbb{S}}_{2}\cup{\mathbb{S}}_{3}=\left\\{\left(k,l\right)\mid{\bm{C}}_{kl}=0\right\\}$
and $\left|{\mathbb{S}}_{2}\cup{\mathbb{S}}_{3}\right|=m+2s$, then when
${\bm{C}}_{\hat{k}\hat{l}}=0$,
${\bm{A}}_{\hat{l}\hat{k}}^{\prime}\left(t\right)\geq\frac{1}{m+2s+{\bm{E}}_{\hat{k}\hat{l}}(t)}$,
where
${\bm{E}}_{\hat{k}\hat{l}}(t)=\sum_{(k,l)\in{\mathbb{S}}_{1}}exp(-{\bm{D}}_{\hat{k},\hat{l},k,l}t)$
and $-{\bm{D}}_{\hat{k},\hat{l},k,l}>0,\forall(k,l)\in{\mathbb{S}}_{1}$ leads
to ${\bm{E}}_{\hat{k}\hat{l}}\left(+\infty\right)=0$. Notice that
$\sum\limits_{k,l}{\bm{A}}_{kl}^{\prime}\left(t\right)=1$, we have

$\frac{1}{m+2s+{\bm{E}}_{\hat{k}\hat{l}}(t)}\leq{\bm{A}}_{\hat{l}\hat{k}}^{\prime}\left(t\right)\leq\frac{1}{m+2s}.$

Therefore,

${\bm{A}}_{kl}^{\prime}\left(+\infty\right)=\left\\{\begin{array}[]{cc}0&,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}\neq{\mathcal{T}}\left({\bm{X}}\right)_{l,:}\\\
\frac{1}{m+2s}&,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}={\mathcal{T}}\left({\bm{X}}\right)_{l,:}\end{array}\right..$

Furthermore,

${\bm{A}}_{kl}\left(+\infty\right)={\bm{A}}_{kl}^{\prime}\left(+\infty\right)+{\bm{A}}_{kl}^{\prime}\left(+\infty\right)=2{\bm{A}}_{kl}^{\prime}\left(+\infty\right)=2{\bm{A}}_{kl}^{\prime}\left(+\infty\right)=\frac{2}{m+2s}=\gamma,$
${\bm{A}}_{kl}\left(+\infty\right)=\left\\{\begin{array}[]{cc}0&,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}\neq{\mathcal{T}}\left({\bm{X}}\right)_{l,:}\\\
\gamma&,{\mathcal{T}}\left({\bm{X}}\right)_{k,:}={\mathcal{T}}\left({\bm{X}}\right)_{l,:}\end{array}\right..$

According to the definition of ${\bm{L}}_{i}$,

${\bm{L}}_{kl}^{*}={\bm{L}}_{kl}(+\infty)=\left\\{\begin{array}[]{cc}0&(k,l)\in{\mathbb{S}}_{1}\\\
\gamma&(k,l)\in{\mathbb{S}}_{2}\\\ -\sum_{l^{\prime}=1,l^{\prime}\neq
k}^{m}{\bm{L}}_{kl^{\prime}}^{*}&k=l\end{array}\right..$

Until now, we have proven that the adaptive regularization part of AIR will
converge at the end of the training, which gives an upper bound of
$\left|{\bm{L}}_{kl}^{*}-{\bm{L}}_{kl}(t)\right|$. Next, we will focus on the
convergence rate of AIR. We separately discuss the rate under three cases in
the formulation above.

If $(k,l)\in{\mathbb{S}}_{2}$ or $k=l$, we denote $D=\min D_{kl}$, according
to the definition of ${\bm{E}}_{kl}(t)$, we have
${\bm{E}}_{kl}(t)\leq\exp(-Dt)$ and

$\displaystyle\left|{\bm{A}}_{kl}^{*}-{\bm{A}}_{kl}(t)\right|$
$\displaystyle=2\left[\frac{1}{m+2s}-\frac{1}{m+2s+{\bm{E}}_{kl}(t)}\right]$
$\displaystyle\leq 2\frac{{\bm{E}}_{kl}(t)}{\left(m+2s\right)^{2}}\leq
2\frac{\frac{m\left(m-1\right)}{2}\exp(-Dt)}{\left(m+2s\right)^{2}}\leq\exp\left(-Dt\right).$

In particular, when $(k,l)\in{\mathbb{S}}_{2}$ we have

$\left|{\bm{L}}_{kl}^{*}-{\bm{L}}_{kl}(t)\right|=\left|{\bm{A}}_{kl}^{*}-{\bm{A}}_{kl}(t)\right|\leq\exp\left(-Dt\right).$

If $(k,l)\in{\mathbb{S}}_{1}$,

$\displaystyle\left|{\bm{L}}_{kl}^{*}-{\bm{L}}_{kl}(t)\right|$
$\displaystyle=\left|{\bm{A}}_{kl}^{*}-{\bm{A}}_{kl}(t)\right|=\left|{\bm{A}}_{kl}(t)\right|$
$\displaystyle\leq\left|\sum_{k^{\prime},l^{\prime}\in{\mathbb{S}}_{1}\cup{\mathbb{S}}_{3}}{\bm{A}}_{k^{\prime}l^{\prime}}(t)\right|$
$\displaystyle=\left|2-\sum_{k^{\prime},l^{\prime}\in{\mathbb{S}}_{2}\cup{\mathbb{S}}_{3}}{\bm{A}}_{k^{\prime}l^{\prime}}(t)\right|$
$\displaystyle=\left|\gamma\left(m+2s\right)-\sum_{k^{\prime},l^{\prime}\in{\mathbb{S}}_{2}\cup{\mathbb{S}}_{3}}{\bm{A}}_{k^{\prime}l^{\prime}}(t)\right|$
$\displaystyle=\left|\sum_{k^{\prime},l^{\prime}\in{\mathbb{S}}_{2}\cup{\mathbb{S}}_{3}}\left(\gamma-{\bm{A}}_{k^{\prime}l^{\prime}}(t)\right)\right|$
$\displaystyle\leq\sum_{k^{\prime},l^{\prime}\in{\mathbb{S}}_{2}\cup{\mathbb{S}}_{3}}\left|\gamma-{\bm{A}}_{k^{\prime}l^{\prime}}(t)\right|$
$\displaystyle=\sum_{k^{\prime},l^{\prime}\in{\mathbb{S}}_{2}\cup{\mathbb{S}}_{3}}\left|{\bm{A}}_{k^{\prime}l^{\prime}}^{*}-{\bm{A}}_{k^{\prime}l^{\prime}}(t)\right|\leq
2\exp(-Dt)/\gamma.$

If $k=l$,

$\left|{\bm{L}}_{kl}^{*}-{\bm{L}}_{kl}(t)\right|=\left|\sum_{l^{\prime}=1,l^{\prime}\neq
k}^{m}\left({\bm{L}}_{kl^{\prime}}^{*}-{\bm{L}}_{kl^{\prime}}(t)\right)\right|\leq
2\left(m-1\right)\exp(-Dt)/\gamma.$

∎

##  Appendix C Deeper and wider are beneficial

We will illustrate that the multi-layer matrix factorization enjoys implicit
low-rank regularization without any special requirement on the initialization
and constraint on shared dimension $r$ (width of the matrix), which overcomes
the sensitivity on initialization in the two-layer matrix factorization [17].
That is, we do not need to estimate $r$ in advance. We design this experiment
to explore the effects of $L$ and $r$ in Figure 10, and we see that the models
perform better with bigger $L$ and $r$. However, bigger $L$ and $r$ also
increase the computation complexity. Thus, we use $L=3$ and
$r=\min\left\\{m,n\right\\}$ as our factorized model.

|
---|---
(a) The Different numbers of the factorized matrix. | (b)Different shared dimensions.
Figure 10: Trajectories of training AIR for inpainting the Barbara image with
randomly missing $50\%$ pixels with (a) different number of factorized matrix,
and (b) different shared common dimension.

## References

  * [1] Joshua B. Tenenbaum, Vin De Silva, and John C. Langford.  A global geometric framework for nonlinear dimensionality reduction.  Science, 290 5500:2319–23, 2000.
  * [2] Emmanuel J. Candès, Michael B. Wakin, and Stephen P. Boyd.  Enhancing sparsity by reweighted $\ell$1 minimization.  Journal of Fourier Analysis and Applications, 14:877–905, 2007\.
  * [3] Weisheng Dong, Guangming Shi, and Xin Li.  Nonlocal image restoration with bilateral variance estimation: A low-rank approach.  IEEE Transactions on Image Processing, 22:700–711, 2013.
  * [4] Remi Lam, Olivier Zahm, Youssef M. Marzouk, and Karen Willcox.  Multifidelity dimension reduction via active subspaces.  SIAM Journal on Scientific Computing, 42:1, 2020.
  * [5] Daniele Bigoni, Youssef M. Marzouk, Clémentine Prieur, and Olivier Zahm.  Nonlinear dimension reduction for surrogate modeling using gradient information.  ArXiv, abs/2102.10351, 2021.
  * [6] Meyer Scetbon, Michael Elad, and Peyman Milanfar.  Deep k-svd denoising.  IEEE Transactions on Image Processing, 30:5944–5955, 2021.
  * [7] Alona Golts, Daniel Freedman, and Michael Elad.  Deep energy: Task driven training of deep neural networks.  IEEE Journal of Selected Topics in Signal Processing, 15:324–338, 2021.
  * [8] Rajaei Khatib, Dror Simon, and Michael Elad.  Learned greedy method (lgm): A novel neural architecture for sparse coding and beyond.  J. Vis. Commun. Image Represent., 77:103095, 2021.
  * [9] Maryam Fazel.  Matrix rank minimization with applications.  PhD thesis, PhD thesis, Stanford University, 2002.
  * [10] Emmanuel J. Candès and Benjamin Recht.  Exact matrix completion via convex optimization.  Foundations of Computational Mathematics, 9:717–772, 2009.
  * [11] Amit Boyarski, Sanketh Vedula, and Alexander M. Bronstein.  Spectral geometric matrix completion.  arXiv: Learning, 2019.
  * [12] Daniel D. Lee and H. Sebastian Seung.  Learning the parts of objects by non-negative matrix factorization.  Nature, 401:788–791, 1999.
  * [13] Andrzej Cichocki, Rafał Zdunek, A. Phan, and Shun-ichi Amari.  Nonnegative matrix and tensor factorizations - applications to exploratory multi-way data analysis and blind source separation.  IEEE Signal Processing Magazine, 25:142–145, 2009.
  * [14] Ajit P. Singh and Geoffrey J. Gordon.  A unified view of matrix factorization models.  In European conference on Machine Learning, 2008.
  * [15] Yu-Xiong Wang and Yu-Jin Zhang.  Nonnegative matrix factorization: A comprehensive review.  IEEE Transactions on Knowledge and Data Engineering, 25:1336–1353, 2013.
  * [16] Sanjeev Arora, Nadav Cohen, W. Hu, and Yuping Luo.  Implicit regularization in deep matrix factorization.  In NeurIPS, 2019.
  * [17] Amy N. Langville, Carl D. Meyer, Russell Albright, James Cox, and David Duling.  Algorithms, initializations, and convergence for the nonnegative matrix factorization.  arXiv: Numerical Analysis, 2014.
  * [18] Suriya Gunasekar, Blake E. Woodworth, Srinadh Bhojanapalli, Behnam Neyshabur, and Nathan Srebro.  Implicit regularization in matrix factorization.  2018 Information Theory and Applications Workshop (ITA), pages 1–10, 2018.
  * [19] Kaiming He, X. Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  2015 IEEE International Conference on Computer Vision (ICCV), pages 1026–1034, 2015.
  * [20] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  In Neural Information Processing Systems, 2018.
  * [21] Simon S. Du, Xiyu Zhai, Barnabás Póczos, and Aarti Singh.  Gradient descent provably optimizes over-parameterized neural networks.  In International Conference on Learning Representations, 2018.
  * [22] Lénaïc Chizat, Edouard Oyallon, and Francis Bach.  On lazy training in differentiable programming.  In Neural Information Processing Systems, 2019.
  * [23] Sanjeev Arora, Simon S. Du, Wei Hu, Zhiyuan Li, Ruslan Salakhutdinov, and Ruosong Wang.  On exact computation with an infinitely wide neural net.  In Neural Information Processing Systems, 2019.
  * [24] Jaehoon Lee, Lechao Xiao, Samuel S. Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jeffrey Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  In Neural Information Processing Systems, 2020.
  * [25] Jiaoyang Huang and Horng-Tzer Yau.  Dynamics of deep neural networks and neural tangent hierarchy.  In International Conference on Machine Learning, 2020.
  * [26] Lénaïc Chizat and Francis Bach.  On the global convergence of gradient descent for over-parameterized models using optimal transport.  In Neural Information Processing Systems, 2018.
  * [27] Grant M. Rotskoff and Eric Vanden-Eijnden.  Parameters as interacting particles: long time convergence and asymptotic error scaling of neural networks.  In Neural Information Processing Systems, 2018.
  * [28] Song Mei, Andrea Montanari, and Phan-Minh Nguyen.  A mean field view of the landscape of two-layer neural networks.  Proceedings of the National Academy of Sciences of the United States of America, 115:201806579, 2018.
  * [29] Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Mean-field theory of two-layers neural networks: dimension-free bounds and kernel limit.  In Conference on Learning Theory, 2019.
  * [30] Lénaïc Chizat and Francis Bach.  Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss.  In Conference on Learning Theory, 2020.
  * [31] Blake Woodworth, Suriya Gunasekar, Jason D. Lee, Edward Moroshko, Pedro Savarese, Itay Golan, Daniel Soudry, and Nathan Srebro.  Kernel and rich regimes in overparametrized models.  In Conference on Learning Theory, 2020.
  * [32] Edward Moroshko, Blake Woodworth, Suriya Gunasekar, Jason D. Lee, Nati Srebro, and Daniel Soudry.  Implicit bias in deep linear classification: Initialization scale vs training accuracy.  In Neural Information Processing Systems, 2020.
  * [33] Arthur Jacot, F. Ged, Franck Gabriel, Berfin cSimcsek, and Clément Hongler.  Deep linear networks dynamics: Low-rank biases induced by initialization scale and l2 regularization.  ArXiv, abs/2106.15933, 2021.
  * [34] Zhemin Li, Zhi-Qin John Xu, Tao Luo, and Hongxia Wang.  A regularised deep matrix factorized model of matrix completion for image restoration.  IET Image Processing, 2022.
  * [35] Chong You, Zhihui Zhu, Qing Qu, and Yi Ma.  Robust recovery via implicit bias of discrepant learning rates for double over-parameterization.  In Neural Information Processing Systems, 2020.
  * [36] L. Rudin, S. Osher, and E. Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: Nonlinear Phenomena, 60:259–268, 1992.
  * [37] Jing Qin, Harlin Lee, Jocelyn T. Chi, Lucas Drumetz, Jocelyn Chanussot, Yifei Lou, and Andrea L. Bertozzi.  Blind hyperspectral unmixing based on graph total variation regularization.  IEEE Transactions on Geoscience and Remote Sensing, 2021.
  * [38] Mujibur Rahman Chowdhury, Jing Qin, and Yifei Lou.  Non-blind and blind deconvolution under poisson noise using fractional-order total variation.  Journal of Mathematical Imaging and Vision, 2020.
  * [39] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Advances in neural information processing systems, pages 5998–6008, 2017.
  * [40] Patrick Mullen, Y. Tong, Pierre Alliez, and Mathieu Desbrun.  Spectral conformal parameterization.  Computer Graphics Forum, 27, 2008.
  * [41] Boris T Polyak.  Some methods of speeding up the convergence of iteration methods.  USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964.
  * [42] Yurii E Nesterov.  A method for solving the convex programming problem with convergence rate o (1/k^ 2).  In Dokl. Akad. Nauk Sssr, volume 269, pages 543–547, 1983.
  * [43] Bao Wang, Tan M Nguyen, Tao Sun, Andrea L Bertozzi, Richard G Baraniuk, and Stanley J Osher.  Scheduled restart momentum for accelerated stochastic gradient descent.  arXiv preprint arXiv:2002.10583, 2020.
  * [44] Bao Wang and Qiang Ye.  Stochastic gradient descent with nonlinear conjugate gradient-style adaptive momentum.  arXiv preprint arXiv:2012.02188, 2020.
  * [45] Tao Sun, Huaming Ling, Zuoqiang Shi, Dongsheng Li, and Bao Wang.  Training deep neural networks with adaptive momentum inspired by the quadratic optimization.  arXiv preprint arXiv:2110.09057, 2021.
  * [46] Tan Nguyen, Richard Baraniuk, Andrea Bertozzi, Stanley Osher, and Bao Wang.  MomentumRNN: Integrating Momentum into Recurrent Neural Networks.  In Advances in Neural Information Processing Systems (NeurIPS 2020), 2020.
  * [47] Hedi Xia, Vai Suliafu, Hangjie Ji, Tan M. Nguyen, Andrea L. Bertozzi, Stanley J. Osher, and Bao Wang.  Heavy ball neural ordinary differential equations.  arXiv preprint arXiv:2010.04840, 2021.
  * [48] Federico Monti, Michael M. Bronstein, and Xavier Bresson.  Geometric matrix completion with recurrent multi-graph neural networks.  In NIPS, 2017.
  * [49] Aanchal Mongia and A. Majumdar.  Drug-target interaction prediction using multi graph regularized nuclear norm minimization.  PLoS ONE, 15, 2020.
  * [50] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In International Conference on Learning Representations, 2015.
  * [51] Jacob Goldberger, Sam T. Roweis, Geoffrey E. Hinton, and Ruslan Salakhutdinov.  Neighbourhood components analysis.  In NIPS, 2004.
  * [52] O. Troyanskaya, M. Cantor, G. Sherlock, P. Brown, T. Hastie, R. Tibshirani, D. Botstein, and R. Altman.  Missing value estimation methods for dna microarrays.  Bioinformatics, 17 6:520–5, 2001.
  * [53] Mingming Yang and Songhua Xu.  A novel patch-based nonlinear matrix completion algorithm for image analysis through convolutional neural network.  Neurocomputing, 389:56–82, 2020.


closed strings, among others. Secondly, since the soft limit can be taken for
any massless particle, it is intriguing to investigate the applicability of
our method in a broader scope, such as scalar effective field theories like
the non-linear sigma model and beyond, as well as amplitudes of theories in
which massless particles couple to fermions.

## Acknowledgments

KZ is supported by NSFC under Grant No. 11805163.

## References

  * (1) H. Elvang and Y.-t. Huang, Scattering Amplitudes, arXiv:1308.1697.
  * (2) C. Cheung, TASI Lectures on Scattering Amplitudes, pp. 571–623.  2018\.  arXiv:1708.03872.
  * (3) R. Britto, F. Cachazo, and B. Feng, New recursion relations for tree amplitudes of gluons, Nucl. Phys. B 715 (2005) 499–522, [hep-th/0412308].
  * (4) R. Britto, F. Cachazo, B. Feng, and E. Witten, Direct proof of tree-level recursion relation in Yang-Mills theory, Phys. Rev. Lett. 94 (2005) 181602, [hep-th/0501052].
  * (5) F. E. Low, Bremsstrahlung of very low-energy quanta in elementary particle collisions, Phys. Rev. 110 (1958) 974–977.
  * (6) S. Weinberg, Infrared photons and gravitons, Phys. Rev. 140 (1965) B516–B524.
  * (7) F. Cachazo and A. Strominger, Evidence for a New Soft Graviton Theorem, arXiv:1404.4091.
  * (8) E. Casali, Soft sub-leading divergences in Yang-Mills amplitudes, JHEP 08 (2014) 077, [arXiv:1404.5551].
  * (9) B. U. W. Schwab and A. Volovich, Subleading Soft Theorem in Arbitrary Dimensions from Scattering Equations, Phys. Rev. Lett. 113 (2014), no. 10 101601, [arXiv:1404.7749].
  * (10) N. Afkhami-Jeddi, Soft Graviton Theorem in Arbitrary Dimensions, arXiv:1405.3533.
  * (11) Z. Bern, S. Davies, and J. Nohle, On Loop Corrections to Subleading Soft Behavior of Gluons and Gravitons, Phys. Rev. D 90 (2014), no. 8 085015, [arXiv:1405.1015].
  * (12) S. He, Y.-t. Huang, and C. Wen, Loop Corrections to Soft Theorems in Gauge Theories and Gravity, JHEP 12 (2014) 115, [arXiv:1405.1410].
  * (13) F. Cachazo and E. Y. Yuan, Are Soft Theorems Renormalized?, arXiv:1405.3413.
  * (14) M. Bianchi, S. He, Y.-t. Huang, and C. Wen, More on Soft Theorems: Trees, Loops and Strings, Phys. Rev. D 92 (2015), no. 6 065022, [arXiv:1406.5155].
  * (15) A. Sen, Subleading Soft Graviton Theorem for Loop Amplitudes, JHEP 11 (2017) 123, [arXiv:1703.00024].
  * (16) D. Nguyen, M. Spradlin, A. Volovich, and C. Wen, The Tree Formula for MHV Graviton Amplitudes, JHEP 07 (2010) 045, [arXiv:0907.2276].
  * (17) C. Boucher-Veronneau and A. J. Larkoski, Constructing Amplitudes from Their Soft Limits, JHEP 09 (2011) 130, [arXiv:1108.5385].
  * (18) L. Rodina, Scattering Amplitudes from Soft Theorems and Infrared Behavior, Phys. Rev. Lett. 122 (2019), no. 7 071601, [arXiv:1807.09738].
  * (19) S. Ma, R. Dong, and Y.-J. Du, Constructing EYM amplitudes by inverse soft limit, JHEP 05 (2023) 196, [arXiv:2211.10047].
  * (20) C. Cheung, K. Kampf, J. Novotny, and J. Trnka, Effective Field Theories from Soft Limits of Scattering Amplitudes, Phys. Rev. Lett. 114 (2015), no. 22 221602, [arXiv:1412.4095].
  * (21) C. Cheung, K. Kampf, J. Novotny, C.-H. Shen, and J. Trnka, On-Shell Recursion Relations for Effective Field Theories, Phys. Rev. Lett. 116 (2016), no. 4 041601.
  * (22) H. Luo and C. Wen, Recursion relations from soft theorems, JHEP 03 (2016) 088, [arXiv:1512.06801].
  * (23) H. Elvang, M. Hadjiantonis, C. R. T. Jones, and S. Paranjape, Soft Bootstrap and Supersymmetry, JHEP 01 (2019) 195, [arXiv:1806.06079].
  * (24) K. Zhou, Tree level amplitudes from soft theorems, JHEP 03 (2023) 021, [arXiv:2212.12892].
  * (25) F.-S. Wei and K. Zhou, Expanding single-trace YMS amplitudes with gauge-invariant coefficients, Eur. Phys. J. C 84 (2024), no. 1 29, [arXiv:2306.14774].
  * (26) C. Hu and K. Zhou, Recursive construction for expansions of tree Yang–Mills amplitudes from soft theorem, Eur. Phys. J. C 84 (2024), no. 3 221, [arXiv:2311.03112].
  * (27) Y.-J. Du and K. Zhou, Multi-trace YMS amplitudes from soft behavior, JHEP 03 (2024) 081, [arXiv:2401.03879].
  * (28) C.-H. Fu, Y.-J. Du, R. Huang, and B. Feng, Expansion of Einstein-Yang-Mills Amplitude, JHEP 09 (2017) 021, [arXiv:1702.08158].
  * (29) B. Feng, X. Li, and K. Zhou, Expansion of Einstein-Yang-Mills theory by differential operators, Phys. Rev. D 100 (2019), no. 12 125012, [arXiv:1904.05997].
  * (30) J. Polchinski, String theory. Vol. 1: An introduction to the bosonic string.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 12, 2007.
  * (31) E. H. Simmons, Dimension-six Gluon Operators as Probes of New Physics, Phys. Lett. B 226 (1989) 132–136.
  * (32) E. H. Simmons, Higher dimension gluon operators and hadronic scattering, Phys. Lett. B 246 (1990) 471–476.
  * (33) P. L. Cho and E. H. Simmons, Looking for gluon substructure at the tevatron, Phys. Lett. B 323 (1994) 401–407, [hep-ph/9307345].
  * (34) H. Kawai, D. C. Lewellen, and S. H. H. Tye, A Relation Between Tree Amplitudes of Closed and Open Strings, Nucl. Phys. B 269 (1986) 1–23.
  * (35) Z. Bern, J. J. M. Carrasco, and H. Johansson, New Relations for Gauge-Theory Amplitudes, Phys. Rev. D 78 (2008) 085011, [arXiv:0805.3993].
  * (36) M. Chiodaroli, M. Günaydin, H. Johansson, and R. Roiban, Scattering amplitudes in $\mathcal{N}=2$ Maxwell-Einstein and Yang-Mills/Einstein supergravity, JHEP 01 (2015) 081, [arXiv:1408.0764].
  * (37) H. Johansson and A. Ochirov, Color-Kinematics Duality for QCD Amplitudes, JHEP 01 (2016) 170, [arXiv:1507.00332].
  * (38) H. Johansson and A. Ochirov, Double copy for massive quantum particles with spin, JHEP 09 (2019) 040, [arXiv:1906.12292].
  * (39) F. Cachazo, S. He, and E. Y. Yuan, Scattering of Massless Particles: Scalars, Gluons and Gravitons, JHEP 07 (2014) 033, [arXiv:1309.0885].
  * (40) M. Chiodaroli, M. Gunaydin, H. Johansson, and R. Roiban, Explicit Formulae for Yang-Mills-Einstein Amplitudes from the Double Copy, JHEP 07 (2017) 002, [arXiv:1703.00421].
  * (41) S. He and Y. Zhang, New Formulas for Amplitudes from Higher-Dimensional Operators, JHEP 02 (2017) 019, [arXiv:1608.08448].
  * (42) Q. Bonnefoy, G. Durieux, and J. Roosmale Nepveu, Higher-derivative relations between scalars and gluons, arXiv:2310.13041.
  * (43) C. Cheung, C.-H. Shen, and C. Wen, Unifying Relations for Scattering Amplitudes, JHEP 02 (2018) 095, [arXiv:1705.03025].
  * (44) K. Zhou and B. Feng, Note on differential operators, CHY integrands, and unifying relations for amplitudes, JHEP 09 (2018) 160, [arXiv:1808.06835].
  * (45) M. Bollmann and L. Ferro, Transmuting CHY formulae, JHEP 01 (2019) 180, [arXiv:1808.07451].
  * (46) Y.-J. Du, B. Feng, and F. Teng, Expansion of All Multitrace Tree Level EYM Amplitudes, JHEP 12 (2017) 038, [arXiv:1708.04514].