# Quadratic Matrix Factorization with Applications to Manifold Learning

Zheng Zhai, Hengchao Chen, and Qiang Sun  Zheng Zhai, Hengchao Chen, and Qiang
Sun are with the Department of Statistical Sciences, University of Toronto
Corresponding author: Qiang Sun. E-mail<EMAIL_ADDRESS>Manuscript
received April 19, 2005; revised August 26, 2015.

###### Abstract

Matrix factorization is a popular framework for modeling low-rank data
matrices. Motivated by manifold learning problems, this paper proposes a
quadratic matrix factorization (QMF) framework to learn the curved manifold on
which the dataset lies. Unlike local linear methods such as the local
principal component analysis, QMF can better exploit the curved structure of
the underlying manifold. Algorithmically, we propose an alternating
minimization algorithm to optimize QMF and establish its theoretical
convergence properties. Moreover, to avoid possible over-fitting, we then
propose a regularized QMF algorithm and discuss how to tune its regularization
parameter. Finally, we elaborate how to apply the regularized QMF to manifold
learning problems. Experiments on a synthetic manifold learning dataset and
two real datasets, including the MNIST handwritten dataset and a cryogenic
electron microscopy dataset, demonstrate the superiority of the proposed
method over its competitors.

###### Index Terms:

Quadratic matrix factorization, alternating minimization, convergence
property, manifold learning.

##  1 Introduction

Matrix factorization has achieved many successes in various applications,
including factor models [1], clustering [2, 3], recommendation system [4],
graph and representation learning [5]. The key idea behind matrix
factorization is that any data matrix admitting a low-rank structure can be
represented as a product of two matrices with smaller dimensions. In a general
form, matrix factorization solves

$\min_{U\in{\mathbb{R}}^{D\times r},{V\in{\mathbb{R}}^{r\times m}}\atop
U,V\in\mathcal{C}}\|X-UV\|_{\rm F}^{2}.$ (1)

where $X\in\mathbb{R}^{D\times m}$ is the data matrix with dimension $D$ and
sample size $m$, $U\in\mathbb{R}^{D\times r}$ and $V\in\mathbb{R}^{r\times m}$
are factors of rank $r$, $\mathcal{C}$ is the feasible set encoding additional
structural information, and $\|\cdot\|_{\rm F}$ is the Frobenius norm.
Building upon (1), many well-known algorithms in the literature can be
obtained by taking specific feasible sets $\mathcal{C}$. For example, if
$\mathcal{C}=\\{V:VV^{T}=I_{r}\\}$ enforces orthonormal constraints on the
rows of $V$, then (1) reduces to principal component analysis (PCA). If
$\mathcal{C}$ enforces non-negative constraints on both $U$ and $V$, or either
$U$ or $V$, then (1) becomes nonnegative matrix factorization (NMF) [6] or
semi-NMF [7] respectively. NMF also shares a strong connection to spectral
clustering [8].

Although matrix factorization (1) has been widely studied in various forms, it
is less explored in the nonlinear case: some rows of $V$ might be nonlinear
functions of other rows of $V$. Nonlinear matrix factorization naturally
arises in manifold learning problems, where data
$\\{x_{i}\\}_{i=1}^{m}\subseteq{\mathbb{R}}^{D}$ are assumed to concentrate
near an unknown $d$-dimensional manifold $\cal M$ and the goal is to recover
$\cal M$ from $\\{x_{i}\\}_{i=1}^{m}$ [9, 10]. Mathematically, assume
$\\{x_{i}\\}_{i=1}^{m}$ are given by

$x_{i}=f(\tau_{i})+\epsilon_{i},\quad i=1,\ldots,m,$

where $f$ is a bijective smooth mapping from an open set ${\cal
T}\subseteq\mathbb{R}^{d}$ to ${\cal M}\subseteq\mathbb{R}^{D}$,
$\tau_{i}\in{\cal T}$ is the $d$-dimensional representation of $x_{i}$, and
$\epsilon_{i}$ is the $i$-th approximation error. Here $f$ and $\tau_{i}$ are
non-identifiable in the sense that
$f(\tau_{i})=f^{\prime}(\tau^{\prime}_{i})$, where $f^{\prime}=f\circ g$ and
$\tau^{\prime}_{i}=g^{-1}(\tau_{i})$ for any diffeomorphism $g$ on $\cal T$,
but $f(\tau_{i})$ is uniquely defined. To find $f(\tau_{i})$, we propose to
minimize the residual sum of squares with respect to $f$ and $\tau_{i}$:

$\min_{f\in{\cal
F},\\{\tau_{i}\\}_{i=1}^{m}\subseteq\mathbb{R}^{d}}\sum_{i=1}^{m}\|x_{i}-f(\tau_{i})\|_{\rm
F}^{2}.$ (2)

Here ${\cal F}=\\{f:\mathbb{R}^{d}\mapsto\mathbb{R}^{D}\\}$ is a prediction
function class. Let $X=(x_{1},\ldots,x_{m})$,
$\Phi=(\tau_{1},\ldots,\tau_{m})$, and
$\Xi(f,\Phi)=(f(\tau_{1}),\ldots,f(\tau_{m}))$. Then (2) can be rewritten as

$\displaystyle\min_{f\in{\cal F},\Phi\in\mathbb{R}^{m\times
d}}\|X-\Xi(f,\Phi)\|_{\rm F}^{2}.$ (3)

If $\cal F$ consists of linear functions only, i.e., ${\cal F}=\\{f\mid
f(\tau)=A\tau,A\in\mathbb{R}^{D\times d}\\}$, then the optimization problem
(3) reduces to the matrix factorization problem (1) with $U=A$, $V=\Phi$, and
$\mathcal{C}=\mathbb{R}^{D\times r}\times\mathbb{R}^{r\times m}$. This is
referred to as linear matrix factorization (LMF). The linear assumption on $f$
can be too restrictive for manifold learning problems because it does not take
the curved structure into account and thus is only applicable to model flat
manifolds, i.e., manifolds that are locally isometric to the Euclidean spaces.
For general manifolds, it is better to consider $f$ as quadratic functions.
Higher-order polynomial functions for $f$ are also possible but tend to
overfit noisy data, rendering poor generalization performance. For the reasons
above, this paper focuses on the quadratic function class

$\begin{gathered}{\cal F}=\\{f(\tau)=c+A\tau+{\cal
B}(\tau,\tau):c\in\mathbb{R}^{D},A\in\mathbb{R}^{D\times d},\\\ \textnormal{
symmetric tensor }\mathcal{B}\in\mathbb{R}^{D\times d\times
d}\\}.\end{gathered}$ (4)

For any quadratic function $f$, there exists a unique matrix
$R\in\mathbb{R}^{D\times(2+3d+d^{2})/2}$ such that $f(\tau)=R\xi(\tau)$, where

$\begin{gathered}\xi(\tau)=[1,\tau^{T},\psi(\tau)^{T}]^{T}\\\
\psi(\tau)=[\tau_{[1]}^{2},\tau_{[1]}\tau_{[2]},...,\tau_{[1]}\tau_{[d]},\tau^{2}_{[2]},...,\tau_{[d]}^{2}]^{T}\in\mathbb{R}^{(d^{2}+d)/2}.\end{gathered}$
(5)

Here $\tau_{[i]}$ denotes the $i$-th coordinate of $\tau$ and $\psi(\cdot)$
maps $\tau$ to a vector consisting of all quadratic and interaction terms of
$\\{\tau_{[i]}\\}_{i=1}^{d}$. Let
$T(\Phi)=(\xi(\tau_{1}),\ldots,\xi(\tau_{m}))$. Then the optimization problem
(3) with the quadratic function class (4) reduces to
$\min_{R,\Phi}\|X-RT(\Phi)\|_{\rm F}^{2}.$ To make $\Phi$ identifiable, we
propose to solve the following optimization problem

$\displaystyle\min_{R,\Phi\atop\Phi\Phi^{T}=I_{d},\Phi{\bf 1}_{m}={\bf
0}}\|X-RT(\Phi)\|_{\rm F}^{2}.$ (6)

This is again a special case of the general matrix factorization problem (1),
and we emphasize that $T(\Phi)$ encodes an implicit constraint that the last
$(d^{2}+d)/2$ rows of $T(\Phi)$ are quadratic functions of its second-
to-$(1+d)$-th rows. Thus, we refer to (6) as quadratic matrix factorization
(QMF).

LMF has been widely studied [6, 11, 12, 13], but the proposed QMF is much more
difficult and has received little attention. In this paper, we propose to
optimize the QMF problem (6) with respect to $\Phi$ and $R$ alternatively.
Specifically, with $\Phi$ fixed, optimizing (6) with respect to $R$ is
equivalent to a linear regression problem. With $R$ fixed, minimizing (6) over
$\Phi$ is a non-convex quadratic projection problem, i.e., projecting a target
point onto the quadratic surface determined by $R$. This quadratic projection
problem can be efficiently solved by an alternating minimization algorithm
when $\|RJ\|_{\rm F}$ is small, where $J=({\bf 0}\
I_{(d^{2}+d)/2})^{T}\in\mathbb{R}^{(2+3d+d^{2})/2\times(d^{2}+d)/2}$. This
motivates us to add a regularizer $\lambda\|RJ\|_{\rm F}^{2}$ to (6) and solve
the regularized QMF problem:

$\displaystyle\min_{R,\Phi\atop\Phi\Phi^{T}=I_{d},\Phi{\bf 1}_{m}={\bf
0}}\ell_{\lambda}(R,\Phi)=\|X-RT(\Phi)\|_{\rm F}^{2}+\lambda\|RJ\|_{\rm
F}^{2}.$ (7)

To solve (7), an alternating minimization algorithm is proposed in Section 4.
We also discuss how to tune $\lambda$ properly.

Our contributions are four-fold. First, motivated by manifold learning
problems, we introduce the quadratic matrix factorization model and propose an
alternating minimization algorithm for solving (6). Second, we establish the
theoretical convergence property of the QMF algorithm. Third, motivated by the
theoretical analysis of the QMF algorithm, we propose a regularized QMF
algorithm and give an adaptive parameter tuning method. Finally, we apply the
regularized QMF algorithm to solve general manifold learning problems.
Numerically, we examine the performance of the proposed method on a synthetic
manifold learning dataset and two real-world datasets, including the MNIST
handwritten dataset and a cryogenic electron microscopy dataset, and
demonstrate the superiority of the proposed method over its competitors.

###  1.1 Related Work and Paper Organization

Our QMF model is different from the problems studied in [14, 15], which are
also referred to as quadratic matrix factorization. They consider
approximating a matrix by a product of multiple low-rank matrices, and some
factor matrix appears twice in the approximation, that is, $X\approx UU^{T}$
or $X\approx AUBU^{T}C$ with $U$ being an unknown factor. They focus on
estimating $U$. In contrast, our paper is motivated by manifold learning
problems and focuses on solving (1) with quadratic constraints: some rows of
$V$ are quadratic functions of the other rows of $V$. Their algorithms cannot
be applied to solve our QMF problem (6).

Many manifold learning methods, such as LLE [16], Isomap [17], Laplacian
eigenmaps [18] and diffusion maps [19], aim to find a lower-dimensional
representation of the dataset while preserving certain geometric structures.
In contrast, our target is to recover the underlying manifold structure on
which the dataset lies. Most algorithms with the same purpose are based on
tangent space estimation [20, 21, 22, 23, 24]. These methods share one common
limitation that they do not take the higher-order smoothness into account.
Compared with the aforementioned methods, the local polynomial approximation
algorithm, which takes higher-order smoothness into account, could achieve a
better convergence rate as shown in [9]. However, [9] and [10] mainly study
the statistical properties of the local polynomial fitting algorithms, leaving
several important computational issues, such as the algorithmic convergence
properties untouched. Our paper addresses these computational issues by
providing algorithmic convergence properties and a regularized QMF algorithm
that potentially avoids over-fitting.

The rest of the paper proceeds as follows. In Section 2, we describe the
alternating minimization algorithm for solving the QMF problem (6). As a key
component, we propose an alternating minimization algorithm to solve the
quadratic projection problem and present its theoretical analysis. We
establish the algorithmic convergence property of QMF in Section 3. In Section
4, we develop a regularized QMF algorithm and discuss the tuning method.
Applications to manifold learning algorithms are given in Section 5, and
numerical experiments are carried out in Section 6. We conclude this paper
with several remarks in Section 7 and leave technical proofs in the Appendix.

###  1.2 Notation

Throughout this paper, we denote by $\|\cdot\|$ and $\|\cdot\|_{\rm F}$ the
spectral norm and the Frobenius norm respectively. For a vector
$v\in\mathbb{R}^{D}$, we use $\|v\|_{1}=\sum_{i=1}^{D}|v_{i}|$ to denote its
$\ell_{1}$ norm. For a matrix $M\in\mathbb{R}^{r\times m}$, denote by
$\sigma_{i}(M)$, $\sigma_{\min}(M)=\sigma_{\min\\{r,m\\}}(M)$, and
$\|M\|_{2,1}=\sum_{i=1}^{r}\|M_{i\cdot}\|$ the $i$-th largest singular value,
the smallest singular value, and the $\ell_{2,1}$ norm of $M$ respectively.
Here $M_{i\cdot}$ is the $i$-th row of $M$. Also, denote by $M^{\dagger}$ the
Moore-Penrose inverse of $M$ and by
$P_{M}=M^{T}(MM^{T})^{\dagger}M\in\mathbb{R}^{m\times m}$ the projection
matrix corresponding to the row subspace of $M$, i.e., the subspace in
$\mathbb{R}^{m}$ spanned by the rows of $M$. Let
$\mathcal{B}\in\mathbb{R}^{D\times d\times d}$ be a tensor, then
$\mathcal{B}(\tau,\eta)\in\mathbb{R}^{D}$ for any
$\tau,\eta\in\mathbb{R}^{d}$. We use $B_{k}\in\mathbb{R}^{d\times d}$ to
denote the $k$-th slice of $\mathcal{B}$, i.e.,
$\mathcal{B}(\tau,\eta)_{k}=\tau^{T}B_{k}\eta=\langle
B_{k},\eta\tau^{T}\rangle$ for all $\tau,\eta\in\mathbb{R}^{d}$ and
$k=1,\ldots,D$. We refer to $\mathcal{B}$ as a symmetric tensor if
$\\{B_{k}\\}_{k=1}^{D}$ are all symmetric. For any $\eta\in\mathbb{R}^{d}$,
define ${\cal B}_{\eta}$ as the action of ${\cal B}$ on the vector $\eta$:

${\cal B}_{\eta}=[B_{1}\eta,\ldots,B_{D}\eta]^{T}\in\mathbb{R}^{D\times d},$
(8)

where $B_{k}$ denotes the $k$-th slice of $\mathcal{B}$. Thus,
$\mathcal{B}(\tau,\eta)={\cal B}_{\eta}\tau={\cal B}_{\tau}\eta$ when
$\mathcal{B}$ is symmetric. Denote by
$\mathcal{B}^{*}(\cdot):\mathbb{R}^{D}\mapsto\mathbb{R}^{d\times d}$ the
adjoint operator of $\mathcal{B}$:

$\mathcal{B}^{*}(c)=\sum_{k=1}^{D}c_{k}B_{k}\in\mathbb{R}^{d\times
d},\quad\forall c\in\mathbb{R}^{D},$

where $B_{k}$ is the $k$-th slice of $\mathcal{B}$. Let ${\bf
1}_{m}=[1,\ldots,1]^{T}\in\mathbb{R}^{m}$ and ${\bf 0}$ be an all-zero matrix,
whose size depends on the context. In addition, let $I_{d}$ be the identity
matrix of size $d\times d$. Given two matrices $A,D\in\mathbb{R}^{d\times d}$,
we use $A\succeq D$ (resp. $A\preceq D$) to indicate that $A-D$ (resp. $D-A$)
is a positive semi-definite matrix.

##  2 Quadratic Matrix Factorization

This section presents an alternating minimization algorithm for solving
quadratic matrix factorization. Given
$\\{x_{i}\\}_{i=1}^{m}\subseteq\mathbb{R}^{D}$, the goal is to solve

$\min_{f\in{\cal
F},\\{\tau_{i}\\}_{i=1}^{m}\subseteq{\mathbb{R}}^{d}}\sum_{i=1}^{m}\|x_{i}-f(\tau_{i})\|^{2}_{\rm
F},$ (9)

where ${\cal F}$ is the quadratic function class (4). Before proceeding to the
algorithm, let us first illustrate the advantages of the quadratic function
class over the linear function class in a swiss roll fitting example. As in
Figure 1, we generate noisy data points near a swiss roll, fit such data using
linear and quadratic functions, and compare the fitted curves with the
underlying truth. It turns out that linear fitting tends to return a polygon,
while quadratic fitting could produce curved lines that recover the underlying
truth better. The difference between linear and quadratic fitting becomes more
significant in the central region, where the curvature is large. This
coincides with the intuition that quadratic fitting performs better because it
takes the curvature into consideration, while linear fitting does not.

Figure 1: A comparison between the linear and quadratic fitting in the swiss
roll fitting problem. The above figures display the fitted curves using linear
and quadratic functions, respectively. The dots represent the raw data and the
dashed lines represent the underlying truth.

To formally present the algorithm, we need some notations. Recall that
$\psi(\cdot)$ and $\xi(\cdot)$ are given by (5). For any symmetric tensor
$\mathcal{B}\in\mathbb{R}^{D\times d\times d}$, there exists a unique matrix
$Q\in\mathbb{R}^{D\times(d^{2}+d)/2}$ such that

$\mathcal{B}(\tau,\tau)=Q\psi(\tau),\quad\forall\tau\in\mathbb{R}^{d}.$ (10)

In particular, the $k$-th row $Q_{k\cdot}$ of $Q$ relates to the $k$-th slice
of $\mathcal{B}$ via the equality $Q_{k\cdot}\psi(\tau)=\tau^{T}B_{k}\tau$. We
shall refer to such $Q$ as the matrix representation of the symmetric tensor
$\mathcal{B}$. Similarly, for any quadratic function $f(\tau)=c+A\tau+{\cal
B}(\tau,\tau)$, there exists a unique matrix
$R\in\mathbb{R}^{D\times(2+3d+d^{2})/2}$ such that
$f(\tau)=R\xi(\tau),\forall\tau\in\mathbb{R}^{d}.$ Indeed, $R=[c,A,Q]$ with
$Q$ being the matrix representation of the symmetric tensor $\cal B$. Let
$X=[x_{1},\ldots,x_{m}]\in\mathbb{R}^{D\times m}$,
$\Phi=[\tau_{1},...,\tau_{m}]\in{\mathbb{R}}^{d\times m}$, and

$\Psi(\Phi)=[\psi(\tau_{1}),...,\psi(\tau_{m})]\in{\mathbb{R}}^{\frac{d^{2}+d}{2}\times
m}.$ (11)

Then the optimization problem (9) is equivalent to

$\displaystyle\min_{R,\Phi}\ell(R,\Phi)=\ell(c,A,Q,\Phi)=\|X-RT(\Phi)\|_{\rm
F}^{2},$ (12)

where

$\displaystyle R=R(c,A,Q)=[c,A,Q],$ (13) $\displaystyle
T(\Phi)=[\xi(\tau_{1}),\ldots,\xi(\tau_{m})]=\left[{\bf
1}_{m},\Phi^{T},\Psi^{T}(\Phi)\right]^{T}.$

This can be viewed as a matrix factorization problem, with constraints given
by (13). The last $(d^{2}+d)/2$ rows of $T(\Phi)$ in the constraints (13) are
determined by the second-to-$(1+d)$th rows of $T(\Phi)$ through the mapping
$\psi$. Recall that $\psi(\tau)$ collects all quadratic and interaction terms
of $\\{\tau_{[i]}\\}_{i=1}^{d}$, thus these constraints specify all possible
quadratic constraints. Problem (12) is thus referred to as QMF. LMF is a
special case of QMF with additional constraints that $Q=0$.

However, (12) suffers from non-identifiability issues. This is because the
minima of (12) are determined by $T(\Phi)$ only through its row space. To see
this, we fix $\Phi$ and consider minimizing $\ell(R,\Phi)$ with respect to $R$
only. The minimizer $\widetilde{R}$ and the product $\widetilde{R}T(\Phi)$ are
given by

$\displaystyle\widetilde{R}=\mathop{\mathrm{argmin}}_{R}\ell(R,\Phi)=XT(\Phi)^{T}(T(\Phi)T(\Phi)^{T})^{\dagger},$
(14)
$\displaystyle\widetilde{R}T(\Phi)=XT(\Phi)^{T}(T(\Phi)T(\Phi)^{T})^{\dagger}T(\Phi)=XP_{T(\Phi)},$

where $M^{\dagger}$ denotes the Moore–Penrose inverse of $M$ and
$P_{T(\Phi)}=T(\Phi)^{T}(T(\Phi)T(\Phi)^{T})^{\dagger}T(\Phi)\in\mathbb{R}^{m\times
m}$. Thus the loss of $\widetilde{R}T(\Phi)$ only depends on the row subspace
of $T(\Phi)$. Substituting (14) into (12), we obtain

$\displaystyle\min_{\Phi}\min_{R}\ell(R,\Phi)=\min_{\Phi}\|X-XP_{T(\Phi)}\|_{\rm
F}^{2}.$

In particular, if $\Phi_{1}$ and $\Phi_{2}$ satisfy
$P_{T(\Phi_{1})}=P_{T(\Phi_{2})}$, then
$\min_{R}\ell(R,\Phi_{1})=\min_{R}\ell(R,\Phi_{2})$ and thus it is impossible
to distinguish $\Phi_{1}$ and $\Phi_{2}$ when optimizing $\eqref{equ:9}$.
Proposition 1 provides concrete transformations on $\Phi$ such that
$\min_{R}\ell(R,\Phi)$ stays the same.

######  Proposition 1.

Suppose $\Phi^{\prime}=Z\Phi+u{\bf 1}_{m}^{T}$ for some invertible matrix
$Z\in\mathbb{R}^{d\times d}$ and some vector $u\in\mathbb{R}^{d}$. Then for
any $R$, there exists $R^{\prime}$ such that
$\ell(R^{\prime},\Phi^{\prime})=\ell(R,\Phi)$. In particular, we have
$\min_{R}\ell(R,\Phi^{\prime})=\min_{R}\ell(R,\Phi)$.

The proof of Proposition 1 is left in the Appendix. To overcome non-
identifiability issues, we add additional constraints $\Phi\Phi^{T}=I_{d}$ and
$\Phi{\bf 1}_{m}=0$ to (12) and solve

$\min_{R,\Phi\atop\Phi\Phi^{T}=I_{d},\Phi{\bf
1}_{m}=0}\ell(R,\Phi)=\|X-RT(\Phi)\|_{\rm F}^{2},$ (15)

where $R=R(c,A,Q)$ and $T(\Phi)$ are defined in (13). Several remarks follow.
First, by Proposition 1, the minima of the constrained problem (15) and the
unconstrained problem (12) are the same. Thus these two problems are
equivalent. Second, by introducing constraints $\Phi\Phi^{T}=I_{d}$, we reduce
the solution space from $\mathbb{R}^{d\times m}$ to the Stiefel manifold,
which potentially helps speed up the optimization procedure [25]. Third, we
still do not have any restrictions on $R$ in (15), so optimizing (15) over $R$
with fixed $\Phi$ is still a regression problem with its solution given by
(14). Fourth, the constraint $\Phi\Phi^{T}=I_{d}$ enforces the scale of $\Phi$
to be neither too large nor too small. Thus, the configuration of the
approximation $RT(\Phi)$ largely depends on $R$ and can be controlled by
regularizing $R$ properly. We will explore this last point in Section 4.

Data: $X=[x_{1},\ldots,x_{m}]\in\mathbb{R}^{D\times m}$

Result: $\Phi=[\tau_{1},\ldots,\tau_{m}]\in\mathbb{R}^{d\times m}$, quadratic
function $f(\tau)=R_{t}\xi(\tau)$.

1

2 initialize $\Phi_{0}\in\mathbb{R}^{d\times m}$ as the top $d$ eigenvectors
of $G=(X-\bar{x}{\bf 1}_{m}^{T})^{T}(X-\bar{x}{\bf 1}_{m}^{T})$;

3 while  _$\|\Phi_{t}^{T}\Phi_{t}-\Phi_{t-1}^{T}\Phi_{t-1}\| >\epsilon$_ do

4 update $R_{t}=\mathop{\mathrm{argmin}}_{R}\ell(R,\Phi_{t-1})$ as in (14);

5 for  _$i=1$ to $m$_ do

6 solve the $i$-th projection problem
$\widetilde{\tau}_{i,t}=\mathop{\mathrm{argmin}}_{\tau\in\mathbb{R}^{d}}\|x_{i}-R_{t}\xi(\tau)\|^{2}$;

7 end for

8 set
$\widetilde{\Phi}_{t}=[\widetilde{\tau}_{1,t},\ldots,\widetilde{\tau}_{m,t}]$;

9 update $\Phi_{t}=Z_{t}\widetilde{\Phi}_{t}(I_{m}-{\bf 1}_{m}{\bf
1}_{m}^{T}/m)$ with $Z_{t}$ given by (17);

10

11 end while

Algorithm 1 An alternating minimization algorithm for quadratic matrix
factorization.

Now we present our first main algorithm. We adopt an alternating minimization
strategy to solve (15). To begin with, we initialize
$\Phi_{0}\in\mathbb{R}^{d\times m}$ as the top $d$ eigenvectors of the gram
matrix $G=(X-\bar{x}{\bf 1}_{m}^{T})^{T}(X-\bar{x}{\bf 1}_{m}^{T})$, where
$\bar{x}=\frac{1}{m}\sum_{i=1}^{m}x_{i}$. During the $t$-th loop, we first fix
$\Phi=\Phi_{t-1}$ and update
$R_{t}=\mathop{\mathrm{argmin}}_{R}\ell(R,\Phi_{t-1})$ as in (14). Next, we
fix $R=R_{t}$ and update
$\widetilde{\Phi}_{t}=\mathop{\mathrm{argmin}}_{\Phi}\ell(R_{t},\Phi)$. This
is a separable problem in the sense that $\widetilde{\Phi}_{t}$ is given by
$\widetilde{\Phi}_{t}=[\widetilde{\tau}_{1,t},\ldots,\widetilde{\tau}_{m,t}]$
with

$\widetilde{\tau}_{i,t}=\mathop{\mathrm{argmin}}_{\tau\in\mathbb{R}^{d}}\|x_{i}-R_{t}\xi(\tau)\|^{2},$
(16)

where $\xi(\tau)$ is defined in (5). We refer to (16) as a quadratic
projection problem because it finds the closest point
$R_{t}\xi(\widetilde{\tau}_{i,t})$ to $x_{i}$ on the quadratic surface
$f_{t}(\tau)=R_{t}\xi(\tau)$. At the end of the $t$-th loop, we set
$\Phi_{t}=Z_{t}\widetilde{\Phi}_{t}(I_{m}-{\bf 1}_{m}{\bf 1}_{m}^{T}/m)$ with
$Z_{t}$ given by

$Z_{t}=(\widetilde{\Phi}_{t}(I_{m}-{\bf 1}_{m}^{T}{\bf
1}_{m}/m)\widetilde{\Phi}_{t}^{T})^{-1/2}.$ (17)

This ensures that the constraints $\Phi_{t}\Phi_{t}^{T}=I_{d}$ and
$\Phi_{t}{\bf 1}_{m}=0$ hold. Given a precision level $\epsilon>0$, we will
repeat the above iterations until the stopping criteria
$\|\Phi_{t}^{T}\Phi_{t}-\Phi_{t-1}^{T}\Phi_{t-1}\|\leq\epsilon$ is met, i.e.,
the row subspace of $\Phi_{t}$ stabilizes. Algorithm 1 summarizes the details.
Till now, the only issue we have not yet addressed is how to solve (16), which
will be discussed in the following subsection.

###  2.1 Quadratic Projection

Since the parameters $R_{t}$ and $x_{i}$ are fixed when solving (16), we omit
the subscripts $i$ and $t$ for simplicity. We rewrite the loss function in
(16) as

$h(\tau)=\|x-c-A\tau-\mathcal{B}(\tau,\tau)\|^{2},$ (18)

where $c,A,{\cal B}$ are determined by $R$ via (10) and (13). Minimizing (18)
with respect to $\tau$ is non-convex, and we consider minimizing the following
surrogate loss instead:

$\displaystyle\min_{\tau,\eta}g(\tau,\eta)=\frac{1}{2}\|x-c-A\tau-{\cal
B}(\tau,\eta)\|^{2}$ $\displaystyle+\frac{1}{2}\|x-c-A\eta-{\cal
B}(\tau,\eta)\|^{2}.$

Note that $g$ is symmetric in $\tau$ and $\eta$, and $g(\tau,\tau)=h(\tau)$.
Starting from $\tau_{0}=\eta_{0}={\bf 0}\in\mathbb{R}^{d}$, we update
$\\{\tau_{s},\eta_{s}\\}$ iteratively in the following manner:

$\left\\{\begin{array}[]{l}\tau_{s}=\mathop{\mathrm{argmin}}_{\tau}g(\tau,\eta_{s-1}),\\\
\eta_{s}=\mathop{\mathrm{argmin}}_{\eta}g(\tau_{s},\eta).\end{array}\right.$
(19)

Upon convergence such that
$(\tau^{*},\eta^{*})=\mathop{\mathrm{argmin}}_{\tau,\eta}g(\tau,\eta)$, we
take $\tau^{*}$ to be the solution to (16). In what follows, we show how (19)
can be efficiently solved and prove that $(\tau_{s},\eta_{s})$ converges to
$(\tau^{*},\tau^{*})$ for some stationary point $\tau^{*}$ of $h$ under
certain conditions.

The update rule (19) can be implemented efficiently as all iterates admit
closed-form solutions. Let us fix $\eta=\eta_{s-1}$ and consider optimizing
$g(\tau,\eta_{s-1})$ over $\tau$. Recall ${\cal B}_{\eta}$ is given by (8) for
any $\eta\in\mathbb{R}^{d}$ and $\mathcal{B}(\tau,\eta)={\cal
B}_{\tau}\eta={\cal B}_{\eta}\tau$. Then $g(\tau,\eta_{s-1})$ can be rewritten
as

$\displaystyle g(\tau,\eta_{s-1})=\frac{1}{2}\|x-c-(A+{\cal
B}_{\eta_{s-1}})\tau\|^{2}$ $\displaystyle+\frac{1}{2}\|x-c-A\eta_{s-1}-{\cal
B}_{\eta_{s-1}}\tau\|^{2}.$

This is a quadratic function of $\tau$, so
$\tau_{s}=\mathop{\mathrm{argmin}}_{\tau}g(\tau,\eta_{s-1})$ admits a closed-
form solution

$\tau_{s}=\Gamma_{\eta_{s-1}}^{-1}\zeta_{\eta_{s-1}},$

where $\Gamma_{\eta_{s-1}}$ and $\zeta_{\eta_{s-1}}$ are

$\displaystyle\Gamma_{\eta_{s-1}}$ $\displaystyle=(A+{\cal
B}_{\eta_{s-1}})^{T}(A+{\cal B}_{\eta_{s-1}})+{\cal B}_{\eta_{s-1}}^{T}{\cal
B}_{\eta_{s-1}},$ (20) $\displaystyle\zeta_{\eta_{s-1}}$
$\displaystyle=(A+{\cal B}_{\eta_{s-1}})^{T}(x-c)+{\cal
B}_{\eta_{s-1}}^{T}(x-c-A\eta_{s-1}).$ (21)

Since $g(\tau,\eta)$ is symmetric in $\tau$ and $\eta$, the dual problem
$\eta_{s}=\mathop{\mathrm{argmin}}_{\eta}g(\tau_{s},\eta)$ can be solved
similarly as

$\eta_{s}=\Gamma_{\tau_{s}}^{-1}\zeta_{\tau_{s}},$

where $\Gamma_{\tau_{s}}$ and $\zeta_{\tau_{s}}$ are given by (20) and (21)
with $\eta_{s-1}$ replaced by $\tau_{s}$.

Now we provide conditions under which $(\tau_{s},\eta_{s})$ converges such
that $\tau^{*}=\lim_{s}\tau_{s}=\lim_{s}\eta_{s}$, and $(\tau^{*},\tau^{*})$
is a first-order stationary point of $g$. This implies that $\tau^{*}$ is a
stationary point of $h$. The key is to analyze the Hessian matrix of $g$:

$H_{g}(\tau,\eta)=\begin{bmatrix}\Gamma_{\eta}&H_{\eta\tau}\\\
H_{\tau\eta}&\Gamma_{\tau}\end{bmatrix},$

where $\Gamma_{\eta}$ and $\Gamma_{\tau}$ are given by (20) with parameter
$\eta$ and $\tau$ respectively and $H_{\eta\tau}$ and $H_{\tau\eta}$ are given
by

$\displaystyle H_{\eta\tau}=H_{\eta\tau}^{T}=$
$\displaystyle-2\mathcal{B}^{*}(x-c-A(\tau+\eta)/2-\mathcal{B}(\tau,\eta))$
(22) $\displaystyle+{\cal B}_{\eta}^{T}(A+{\cal B}_{\tau})+(A+{\cal
B}_{\eta})^{T}{\cal B}_{\tau}.$

Here $B_{\eta}$ and $B_{\tau}$ are defined in (8), and
$\mathcal{B}^{*}(\cdot)$ represents the adjoint operator of $\mathcal{B}$.
$\Gamma_{\eta}$ is always positive definite when $\sigma_{d}(A)>0$ because

$\Gamma_{\eta}=\frac{1}{2}A^{T}A+(\frac{1}{\sqrt{2}}A+\sqrt{2}{\cal
B}_{\eta})^{T}(\frac{1}{\sqrt{2}}A+\sqrt{2}{\cal
B}_{\eta})\succeq\frac{1}{2}A^{T}A.$ (23)

The following theorem shows that if $(\tau_{s},\eta_{s})$ is bounded, then
under certain conditions on $\mathcal{B}$, we have $(\tau_{s},\eta_{s})$
converges, $\tau^{*}=\lim_{s\to\infty}\tau_{s}=\lim_{s\to\infty}\eta_{s}$, and
$(\tau^{*},\tau^{*})$ is a stationary point of $g$.

######  Theorem 2.

Suppose $\sigma_{d}(A)>0$ and define
$\mathcal{S}_{\alpha}=\\{\tau\in\mathbb{R}^{d}\mid\|\tau\|\leq\alpha\\}$ for
some $\alpha>0$. Denote by $B_{k}$ the $k$-th slice of $\mathcal{B}$ for
$k=1,\ldots,D$. If $\mathfrak{b}=\max_{k}\sigma_{1}(B_{k})$ satisfies

$\displaystyle(2\|x-c\|_{1}+4\alpha\|A\|_{2,1})\mathfrak{b}+3D\alpha^{2}\mathfrak{b}^{2}\leq\sigma^{2}_{d}(A)/4,$
(24)

then $H_{g}(\tau,\eta)$ is positive definite and
$\sigma_{\min}(H_{g}(\tau,\eta))\geq\sigma_{d}^{2}(A)/4$ for all
$\tau,\eta\in\mathcal{S}_{\alpha}$. If the sequence $(\tau_{s},\eta_{s})$
obtained by (19) falls into the region
$\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$ for sufficiently large $s$,
then $(\tau_{s},\eta_{s})$ converges,
$\tau^{*}=\lim_{s\to\infty}\tau_{s}=\lim_{s\to\infty}\eta_{s}$, and
$(\tau^{*},\tau^{*})$ is a stationary point of $g$.

######  Proof.

Since $\sigma_{d}(A)>0$, by (23), $\Gamma_{\eta}$ is positive definite with
$\sigma_{\min}(\Gamma_{\eta})\geq\sigma_{d}^{2}(A)/2$. By symmetry, such
property holds for $\Gamma_{\tau}$ as well. As for $H_{\eta\tau}$ and
$H_{\tau\eta}$, we can use (24) to show that

$\sigma_{1}(H_{\eta\tau})=\sigma_{1}(H_{\tau\eta})\leq\sigma_{d}^{2}(A)/4,$
(25)

holds for any $\tau,\eta\in\mathcal{S}_{\alpha}$. In particular, for any
$\tau,\eta\in{\cal S}_{\alpha}$, we have

$\begin{array}[]{ll}\sigma_{1}({\cal B}^{*}({\cal B}(\tau,\eta)))\leq
D\alpha^{2}\mathfrak{b}^{2},&\sigma_{1}({\cal
B}_{\eta}^{T}A)\leq\alpha\|A\|_{2,1}\mathfrak{b},\\\ \sigma_{1}(A^{T}{\cal
B}_{\tau})\leq\alpha\|A\|_{2,1}\mathfrak{b},&\sigma_{1}({\cal
B}_{\eta}^{T}{\cal B}_{\tau})\leq D\alpha^{2}\mathfrak{b}^{2},\end{array}$

Similarly,

$\displaystyle\sigma_{1}({\cal B}^{*}(x-c-A(\tau+\eta)/2))$
$\displaystyle\leq$
$\displaystyle(\|x-c\|_{1}+\|A(\tau+\eta)/2\|_{1})\mathfrak{b}\leq(\|x-c\|_{1}+\alpha\|A\|_{2,1})\mathfrak{b}.$

Combining these inequalities with (22), we obtain

$\displaystyle\begin{aligned}
\sigma_{1}(H_{\eta\tau})&=\sigma_{1}(H_{\tau\eta})\\\
&\leq(2\|x-c\|_{1}+4\alpha\|A\|_{2,1})\mathfrak{b}+3D\alpha^{2}\mathfrak{b}^{2}.\end{aligned}$

Then (25) follows from the condition (24). To proceed, we decompose
$H_{g}(\tau,\eta)$ into the following summation of a positive definite matrix
and a symmetric matrix:

$H_{g}(\tau,\eta)=\begin{bmatrix}\Gamma_{\eta}&{\bf 0}\\\ {\bf
0}&\Gamma_{\tau}\end{bmatrix}+\begin{bmatrix}{\bf 0}&H_{\eta\tau}\\\
H_{\tau\eta}&{\bf 0}\end{bmatrix}.$ (26)

For the second term, it follows from (25) that

$\sigma_{1}\left(\begin{bmatrix}{\bf 0}&H_{\eta\tau}\\\ H_{\tau\eta}&{\bf
0}\end{bmatrix}\right)=\sigma_{1}(H_{\eta\tau})\leq\sigma_{d}^{2}(A)/4,$

for any $\tau,\eta\in\mathcal{S}_{\alpha}$. Since
$\min\\{\sigma_{\min}(\Gamma_{\tau}),\sigma_{\min}(\Gamma_{\eta})\\}\geq\sigma_{d}^{2}(A)/2$,
the first term in (26) is positive definite, $H_{g}(\tau,\eta)$ is positive
definite, and

$\displaystyle\sigma_{\min}(H_{g}(\tau,\eta))\geq$
$\displaystyle\min\\{\sigma_{\min}(\Gamma_{\tau}),\sigma_{\min}(\Gamma_{\eta}\\}-\sigma_{1}(H_{\eta\tau})$
(27) $\displaystyle\geq$
$\displaystyle\sigma_{d}^{2}(A)/4,\qquad\forall\tau,\eta\in\mathcal{S}_{\alpha}.$

Since the region $\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$ is convex,
(27) implies the strong convexity of $g$ over
$\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$.

Suppose that the sequence $(\tau_{s},\eta_{s})$ generated by (19) falls into
the region $\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$ for sufficiently
large $s$. Since $g$ is smooth and strongly convex over the region
$\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$, it is well-known that
$(\tau_{s},\eta_{s})$ would converge to the unique first-order stationary
point of $g$ in $\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$ [26]. Next,
we will prove $\lim_{s\to\infty}\tau_{s}=\lim_{s\to\infty}\eta_{s}$ by
contradiction. In specific, define $\tau^{*}=\lim_{s\to\infty}\tau_{s}$ and
$\eta^{*}=\lim_{s\to\infty}\eta_{s}$ and assume $\tau^{*}\neq\eta^{*}$. Since
$(\tau^{*},\eta^{*})$ is a first-order stationary point of $g$ in
$\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$, by symmetry, we know
$(\eta^{*},\tau^{*})$ is also a first-order stationary point of $g$ in
$\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$. However, this contradicts
the uniqueness of the first-order stationary point of $g$ in
$\mathcal{S}_{\alpha}\times\mathcal{S}_{\alpha}$.

Figure 2:  An illustration of the effect of $\mathfrak{b}$ on the convergence
behaviors of $\\{(\tau_{s},\eta_{s})\\}_{s=1}^{\infty}$ generated by (19). The
left column corresponds to the setting $(x,c,A,20{\cal B})$ and the right
column corresponds to the setting $(x,c,A,30{\cal B})$, where the parameters
are set as follows: $A=[-0.8979,1.0086,-0.5422]^{T},{\cal
B}=[0.7817,-1.4908,-0.3679],c=[0.4171,0.9176,0.1759]^{T}$, and
$x=[0.2561,0.7500,0.0099]^{T}$. The first row displays $g(\tau_{s},\eta_{s})$
with the background surface representing the function $g$. The dotted line
stands for the case $\tau=\eta$. The second row displays $h(\tau_{s})$ under
the two settings.

Suppose $(\tau_{s},\eta_{s})\subseteq{\cal S}_{\alpha}\times{\cal S}_{\alpha}$
for sufficiently large $s$. Theorem 2 proves the convergence of
$(\tau_{s},\eta_{s})$ by establishing the strong convexity of $g$ over ${\cal
S}_{\alpha}\times{\cal S}_{\alpha}$ under condition (24). Condition (24) holds
when both $\alpha$ and $\|x-c\|_{1}$ are small or
$\mathfrak{b}=\max_{k}\sigma_{1}(B_{k})$ is small. In particular, condition
(24) holds if $\mathfrak{b}\leq\mathfrak{b}_{0}$, where

$\displaystyle\mathfrak{b}_{0}=$
$\displaystyle\frac{-\|x-c\|_{1}-2\alpha\|A\|_{2,1}}{3D\alpha^{2}}$
$\displaystyle+\frac{\sqrt{(\|x-c\|_{1}+2\alpha\|A\|_{2,1})^{2}+3D\alpha^{2}\sigma_{d}^{2}(A)/4}}{3D\alpha^{2}}.$

It trivially holds when $\mathfrak{b}=0$, which corresponds to LMF. In
general, by the relationship (10) between $\mathcal{B}$ and $Q$, we have
$\sigma_{1}(B_{k})\leq\|B_{k}\|_{\rm F}\leq\|Q_{k\cdot}\|$, where $Q_{k\cdot}$
is the $k$-th row of $Q$. Thus, $\mathfrak{b}$ is small if $\|Q\|_{\rm F}$ is
small.

To conclude this subsection, we use an example with $D=3,d=1$ to illustrate
the effect of $\mathfrak{b}$ on the convergence of
$\\{(\tau_{s},\eta_{s})\\}_{s=1}^{\infty}$. We consider two different settings
of $g$ and display the sequence $\\{(\tau_{s},\eta_{s})\\}_{s=1}^{\infty}$ in
Figure 2. The left panel has a smaller $\mathfrak{b}$ than the right panel.
The first row of Figure 2 shows that $(\tau_{s},\eta_{s})$ converges with
$\lim\tau_{s}=\lim\eta_{s}$ in the left panel while $(\tau_{s},\eta_{s})$
converges with $\lim\tau_{s}\neq\lim\eta_{s}$ in the right panel. In addition,
the second row shows that $\tau_{s}$ in the left panel converges to the global
minimum of $h$ while $\tau_{s}$ in the right panel does not. The differences
between these two cases lie in the different values of $\mathfrak{b}$. In
particular, if $\mathfrak{b}$ is too large as in the right panel,
$(\tau_{s},\eta_{s})$ would not converge to the same limiting point, that is
$\lim\tau_{s}\neq\lim\eta_{s}$. In this case, $\eta_{s}$ is not guaranteed to
converge to the global minimum of $h$ in general.

##  3 Convergence Properties

This section establishes the convergence property for Algorithm 1. In
particular, we provide conditions under which a limit point of the update
sequence is a stationary point of (12).

######  Theorem 3.

Denote by $\\{(R_{t},\widetilde{\Phi}_{t},\Phi_{t})\\}_{t=1}^{\infty}$ the
update sequence generated by Algorithm 1. For any sub-sequence
$\\{t_{j}\\}_{j=1}^{\infty}$ such that $\lim_{j}\Phi_{t_{j}}=\Phi^{*}$,
$\lim_{j}\Phi_{t_{j}-1}=\Phi^{**}$,
$\lim_{j}\widetilde{\Phi}_{t_{j}}=\widetilde{\Phi}^{**}$, if
$\sigma_{\min}(T(\Phi^{*})T(\Phi^{*})^{T})>0$ and
$\sigma_{\min}(T(\Phi^{**})T(\Phi^{**})^{T})>0$, then $R_{t_{j}+1}$ and
$R_{t_{j}}$ also converge. Denote by $R^{*}=\lim_{j}R_{t_{j}+1}$ and
$R^{**}=\lim_{j}R_{t_{j}}$. Also, we have
$\widetilde{\Phi}^{**}=\mathop{\mathrm{argmin}}_{\Phi}\ell(R^{**},\Phi)$. If
we assume

$\ell(R^{**},\Phi^{**})-\ell(R^{**},\widetilde{\Phi}^{**})\geq\gamma\|\Phi^{**}-\widetilde{\Phi}^{**}\|_{\rm
F}^{2}$ (28)

holds for some constant $\gamma>0$, then
$\Phi^{**}=\widetilde{\Phi}^{**}=\Phi^{*}$ and $R^{**}=R^{*}$ and the
accumulation point $(R^{*},\Phi^{*})$ satisfies the first order Karush-Kuhn-
Tucker (KKT) condition of the problem (12).

######  Proof.

By Algorithm 1, the update sequence
$\\{(R_{t},\widetilde{\Phi}_{t},\Phi_{t})\\}_{t=1}^{\infty}$ satisfies

$\begin{cases}R_{t}=\mathop{\mathrm{argmin}}_{R}\ell(R,\Phi_{t-1}),\\\
\widetilde{\Phi}_{t}=\mathop{\mathrm{argmin}}_{\Phi}\ell(R_{t},\Phi),\\\
\Phi_{t}=Z_{t}\widetilde{\Phi}_{t}(I_{m}-{\bf 1}_{m}{\bf
1}_{m}^{T}/m),\end{cases}$

where $Z_{t}$ is defined by (17). It implies that $\ell(R_{t+1},\Phi_{t})$ is
monotone decreasing:

$\displaystyle\ell(R_{t+1},{\Phi}_{t})=$
$\displaystyle\min_{R}\ell(R,\Phi_{t})=\min_{R}\ell(R,\widetilde{\Phi}_{t})$
(29) $\displaystyle\leq$
$\displaystyle\ell(R_{t},\widetilde{\Phi}_{t})\leq\ell(R_{t},\Phi_{t-1}),$

where the second equality follows from Proposition 1. By the monotone
convergence theorem and the fact that $\ell(R_{t+1},\Phi_{t})\geq 0$, both
$\ell(R_{t},\Phi_{t-1})$ and $\ell(R_{t},\widetilde{\Phi}_{t})$ converge to
the same value.

Consider the sub-sequence $\\{t_{j}\\}_{j=1}^{\infty}$ such that
$\lim_{j}\Phi_{t_{j}}=\Phi^{*}$ and
$\sigma_{\min}(T(\Phi^{*})T(\Phi^{*})^{T})=\gamma_{1}>0$. Then
$\sigma_{\min}(T(\Phi)T(\Phi)^{T})\geq\gamma_{1}/2>0$ for any
$\Phi\in\mathcal{S}^{*}_{\epsilon_{0}}=\\{\Phi\mid\|\Phi-\Phi^{*}\|_{\rm
F}\leq\epsilon_{0}\\}$ with some sufficiently small $\epsilon_{0}$. Thus,
$XT(\Phi)^{T}(T(\Phi)T(\Phi)^{T})^{-1}$ is well-defined for
$\Phi\in\mathcal{S}^{*}_{\epsilon_{0}}$ and is a continuous function of $\Phi$
over $\mathcal{S}^{*}_{\epsilon_{0}}$. Since $\lim_{j}\Phi_{t_{j}}=\Phi^{*}$,
we know $\Phi_{t_{j}}\in\mathcal{S}^{*}_{\epsilon_{0}}$ for sufficiently large
$j$. By continuity, we have

$\displaystyle\lim_{j}R_{t_{j}+1}=$
$\displaystyle\lim_{j}XT(\Phi_{t_{j}})^{T}(T(\Phi_{t_{j}})T(\Phi_{t_{j}})^{T})^{\dagger}$
$\displaystyle=$ $\displaystyle
XT(\Phi^{*})^{T}(T(\Phi^{*})T(\Phi^{*})^{T})^{-1}.$

where we use the definition (14) of $R_{t_{j}+1}$. Denote by
$R^{*}=\lim_{j}R_{t_{j}+1}$.

Let us further assume that the sub-sequence $\\{t_{j}\\}_{j=1}^{\infty}$
satisfies $\lim_{j}\Phi_{t_{j}-1}=\Phi^{**}$,
$\lim_{j}\widetilde{\Phi}_{t_{j}}=\widetilde{\Phi}^{**}$, and
$\sigma_{\min}(T(\Phi^{**})T(\Phi^{**})^{T})>0$. Then $R_{t_{j}}$ converges by
the same argument and denote by $R^{**}=\lim_{j}R_{t_{j}}$. Moreover, we claim
that $\widetilde{\Phi}^{**}=\mathop{\mathrm{argmin}}_{\Phi}\ell(R^{**},\Phi)$.
Otherwise, there exists $\Phi^{\prime}$ such that
$\ell(R^{**},\Phi^{\prime})<\ell(R^{**},\widetilde{\Phi}^{**})$. By
continuity, there exists some $j_{0}$ such that
$\ell(R_{t_{j_{0}}},\Phi^{\prime})<\ell(R^{**},\widetilde{\Phi}^{**})$, which
contradicts to the relationship (29). Furthermore, we assume (28) holds for
some constant $\gamma>0$.

Now we are in a position to prove that
$\Phi^{**}=\widetilde{\Phi}^{**}=\Phi^{*}$, $R^{**}=R^{*}$, and
$(R^{*},\Phi^{*})$ satisfies the KKT condition. By continuity, we have

$\displaystyle\ell(R^{**},\Phi^{**})=$
$\displaystyle\lim_{j}\ell(R_{t_{j}},\Phi_{t_{j}-1})$ $\displaystyle=$
$\displaystyle\lim_{j}\ell(R_{t_{j}},\widetilde{\Phi}_{t_{j}})=\ell(R^{**},\widetilde{\Phi}^{**}),$

where the second equality follows from (29). By (28), we have
$\Phi^{**}=\widetilde{\Phi}^{**}$, or equivalently
$\|\Phi_{t_{j}-1}-\widetilde{\Phi}_{t_{j}}\|_{\rm F}^{2}\to 0$ as
$j\to\infty$. Since
$\Phi_{t_{j}}=\mathop{\mathrm{argmin}}_{\Phi\Phi^{T}=I_{d},\Phi{\bf
1}_{m}=0}\|\Phi-\widetilde{\Phi}_{t_{j}}\|_{\rm F}^{2}$, we have
$\|\Phi_{t_{j}}-\widetilde{\Phi}_{t_{j}}\|_{\rm
F}^{2}\leq\|\Phi_{t_{j}-1}-\widetilde{\Phi}_{t_{j}}\|_{\rm F}^{2}$ and thus
$\|\Phi_{t_{j}}-\widetilde{\Phi}_{t_{j}}\|_{\rm F}^{2}\to 0$ as $l\to\infty$.
Then the fact $\Phi_{t_{j}}\to\Phi^{*}$ implies that
$\Phi^{**}=\widetilde{\Phi}^{**}=\Phi^{*}$. By (Proof) and its $R^{**}$
version, we have $R^{**}=R^{*}$. Finally, by taking the limit on the following
optimality conditions:

$\begin{cases}\nabla_{\Phi}\ell(R_{t_{j}},\Phi)|_{\Phi=\widetilde{\Phi}_{t_{j}}}=0,\\\
\Phi_{t_{j}}=Z_{t_{j}}\widetilde{\Phi}_{t_{j}}(I_{m}-{\bf 1}_{m}{\bf
1}_{m}^{T}/m),\\\
\nabla_{R}\ell(R,\Phi_{t_{j}})|_{R=R_{t_{j}+1}}=0,\end{cases}$

we prove that $\\{R^{*},\Phi^{*}\\}$ satisfies the KKT conditions of the
problem (12).

Let us discuss the assumptions in Theorem 3. First, the feasible set
$\mathcal{G}=\\{\Phi\in\mathbb{R}^{d\times m}\mid\Phi\Phi^{T}=I_{d},\Phi{\bf
1}_{m}=0\\}$ is compact, so there always exists a sub-sequence
$\\{t_{j}\\}_{j=1}^{\infty}$ such that $\\{\Phi_{t_{j}}\\}_{j=1}^{\infty}$ and
$\\{\Phi_{t_{j}-1}\\}_{j=1}^{\infty}$ converge. If we assume
$\widetilde{\Phi}_{t}$ is bounded, then we can choose
$\\{t_{j}\\}_{j=1}^{\infty}$ such that $\widetilde{\Phi}_{t_{j}}$ also
converges. Second, the assumption
$\sigma_{\min}(T(\Phi^{*})T(\Phi^{*})^{T})>0$ and
$\sigma_{\min}(T(\Phi^{**})T(\Phi^{**})^{T})>0$ are mild, especially when $m$
is much larger than $d$. Finally, the assumption (28) characterizes the
$\gamma$-strong convexity of $\ell(R^{**},\Phi)$ with respect to $\Phi$
surrounding the minimizer $\widetilde{\Phi}^{**}$. Theorem 2 provides
conditions under which $g(\tau,\eta)$ is strongly convex over a bounded set,
which implies the strong convexity of $\ell(R^{**},\Phi)$ with respect to
$\Phi$ over a bounded set.

##  4 Regularized Quadratic Matrix Factorization

Motivated by the discussion at the end of Section 2.1, we propose a
regularized quadratic matrix factorization (RQMF) method. Specifically, we add
a regularizer $\lambda\|Q\|_{\rm F}$ to the original problem (15) to prevent
$\|Q\|_{\rm F}$ from being too large. This leads to the following RQMF
problem:

$\min_{R,\Phi\atop\Phi\Phi^{T}=I_{d},\Phi{\bf
1}_{m}=0}\ell_{\lambda}(R,\Phi)=\|X-RT(\Phi)\|^{2}_{\rm F}+\lambda\|RJ\|_{\rm
F}^{2},$ (30)

where $R$ and $T(\Phi)$ are given by (13),

$J=\begin{bmatrix}0&I_{(d^{2}+d)/2}\end{bmatrix}^{T}\in\mathbb{R}^{\frac{2+3d+d^{2}}{2}\times\frac{d^{2}+d}{2}},\textnormal{
and }RJ=Q.$ (31)

Since the columns of $X$ are assigned equal weights in (30), we denote it as
RQMF-E in the experiment section to distinguish it from its kernel version. In
the following, we can still denote it as RQMF when not causing confusion.

As shall be seen later, RQMF with a proper $\lambda$ avoids over-fitting and
thus enjoys a better generalization performance. In the rest of this section,
we will first provide an alternating minimization algorithm to solve (30) and
then elaborate how to tune $\lambda$ properly. Also, we show that no matter
how $\lambda$ is chosen, RQMF always outperforms LMF in terms of memorization
properties.

Data: $X=[x_{1},\ldots,x_{m}]\in\mathbb{R}^{D\times m}$

Result: $\Phi=[\tau_{1},\ldots,\tau_{m}]\in\mathbb{R}^{d\times m}$, quadratic
function $f(\tau)=R_{t}\xi(\tau)$.

1 initialize $\Phi_{0}\in\mathbb{R}^{d\times m}$ as the top $d$ eigenvectors
of $G=(X-\bar{x}{\bf 1}_{m}^{T})^{T}(X-\bar{x}{\bf 1}_{m}^{T})$;

2 while  _$\|\Phi_{t}^{T}\Phi_{t}-\Phi_{t-1}^{T}\Phi_{t-1}\| >\epsilon$_ do

3 update $R_{t}=\mathop{\mathrm{argmin}}_{R}\ell_{\lambda}(R,\Phi_{t-1})$ as
in (32);

4 for  _$i=1$ to $m$_ do

5 solve the $i$-th projection problem
$\widetilde{\tau}_{i,t}=\mathop{\mathrm{argmin}}_{\tau\in\mathbb{R}^{d}}\|x_{i}-R_{t}\xi(\tau)\|^{2}$;

6 end for

7 set
$\widetilde{\Phi}_{t}=[\widetilde{\tau}_{1,t},\ldots,\widetilde{\tau}_{m,t}]$;

8 update $\Phi_{t}=Z_{t}\widetilde{\Phi}_{t}(I_{m}-{\bf 1}_{m}{\bf
1}_{m}^{T}/m)$ with $Z_{t}$ given by (17);

9

10 end while

Algorithm 2 Regularized Quadratic Matrix Factorization

To solve (30), we adopt the same alternating minimization strategy described
in Algorithm 1. When $R$ is fixed, minimizing $\ell_{\lambda}(R,\Phi)$ with
respect to $\Phi$ is equivalent to minimizing $\ell(R,\Phi)$ with respect to
$\Phi$, since the regularizer $\lambda\|RJ\|_{\rm F}^{2}$ is independent of
$\Phi$. It thus reduces to the quadratic projection problem discussed in
Section 2.1. On the other hand, when $\Phi$ is fixed, minimizing
$\ell_{\lambda}(R,\Phi)$ with respect to $R$ is a ridge regression problem,
and the solution can be given in closed form as

$\displaystyle\widetilde{R}=$
$\displaystyle\mathop{\mathrm{argmin}}_{R}\ell_{\lambda}(R,\Phi)$ (32)
$\displaystyle=$ $\displaystyle XT(\Phi)^{T}(T(\Phi)T(\Phi)^{T}+\lambda
JJ^{T})^{-1}.$

Here we use the observation that $T(\Phi)T(\Phi)^{T}+\lambda JJ^{T}$ is
invertible as shown in Lemma 4 below. Therefore, to solve (30), it suffices to
replace (14) in Algorithm 1 by (32), which gives us Algorithm 2, the RQMF
algorithm.

######  Lemma 4.

Suppose $\Phi\Phi^{T}=I_{d}$ and $\Phi{\bf 1}_{m}=0$. If $\lambda>0$, then
$T(\Phi)T(\Phi)^{T}+\lambda JJ^{T}$ is positive definite. Furthermore,
$J^{T}(T(\Phi)T(\Phi)^{T}+\lambda JJ^{T})^{-1}J$ is also positive definite.

The proof of Lemma 4 is left in the Appendix. The following proposition shows
that no matter how $\lambda$ is chosen, RQMF memorizes the data better than
LMF. Recall that LMF corresponds to RQMF with $\lambda\to\infty$, or
equivalently $Q=RJ=0$.

######  Proposition 5.

For any $\lambda>0$, the RQMF that solves (30) memorizes the data better than
LMF in the following sense:

$\|X-R^{*}T(\Phi^{*})\|_{\rm F}^{2}\leq\|X-R^{\prime}T(\Phi^{\prime})\|_{\rm
F}^{2},$ (33)

where
$(R^{\prime},\Phi^{\prime})=\mathop{\mathrm{argmin}}_{R\in\Omega,\Phi}\ell_{\lambda}(R,\Phi)$
with $\Omega=\\{R\mid R=[c,A,Q],Q=0\\}$ is the solution of LMF and
$(R^{*},\Phi^{*})=\mathop{\mathrm{argmin}}_{R,\Phi}\ell_{\lambda}(R,\Phi)$ is
the solution of RQMF.

The proof of Proposition 5 is left in the Appendix.

###  4.1 Tuning Parameter Selection

Figure 3: Illustration of the effect of $\lambda$ for fitting a circle. This
figure displays the generated data and the locally fitted curves with
$\lambda=0.1,0.01,0$ (from left to right). In the last three figures, the
rhombuses stand for the target samples, and the asterisks represent the place
where the target samples are projected to fitted curves.

This section discusses how to tune $\lambda$. Before presenting our new tuning
method, let us first illustrate the effect of $\lambda$ in Figure 3. We
generate 240 data points uniformly on the unit circle and then manually add
normal noises obeying ${\cal N}(0,0.1^{2}I)$. For each target sample, we fit a
curve around the target data using the nearest 40 data points and then project
the target data onto the fitted curve. To fit the curve, we use the RQMF
algorithm with $\lambda=0.1$, 0.01, and 0. Figure 3 shows that the RQMF
algorithm with $\lambda=0.01$ achieves the best performance. When
$\lambda=0.1$ is too large, the RQMF algorithm behaves like linear matrix
factorization and tends to use straight lines as the fitted curves. When
$\lambda=0$ is too small, the RQMF algorithm tends to overfit data with
excessively curved lines. Therefore, it is important to pick a proper
$\lambda$.

In what follows, we will describe a new adaptive tuning method. Recall that
when $\Phi$ is fixed, $\widetilde{R}$ in (32) is a function of $\lambda$.
Define $s(\lambda)=\|\widetilde{R}(\lambda)J\|_{\rm F}^{2}$ and denote by
$s^{\prime}(\lambda)$ and $s^{\prime\prime}(\lambda)$ the corresponding first
and second derivatives.

Figure 4: The fitting error of RQMF against different $\lambda$,
$\delta(\lambda)=-s^{\prime}(\lambda)$, and $s(\lambda)$. Here we take
$\lambda\in[0,0.1]$ and compute $\delta(\lambda)$ and $s(\lambda)$. To compare
these three curves in the same horizontal axis $[0,0.1]$, we shift and re-
scale $\delta(\lambda)$ and $s(\lambda)$ via the function:
$f(x)=0.1\cdot\frac{x-\min(x)}{\max(x)-\min(x)}$.

Proposition 6 shows that $s^{\prime}(\lambda)<0$ and
$s^{\prime\prime}(\lambda)>0$ when $s(\lambda)>0$ and $\lambda>0$. This
implies that $s(\lambda)$ is a decreasing function of $\lambda$ while
$s^{\prime}(\lambda)$ is a strictly increasing function of $\lambda$. In
particular, the root $\lambda=s^{\prime-1}(-\delta)$ is unique and can be
easily found via the bisection method. We propose to pick
$\lambda=s^{\prime-1}(-\delta)$ for a prescribed $\delta>0$.

The proposed tuning method is reasonable in the following sense. Recall that
$s(\lambda)$ is the quantity that the regularizer $\|RJ\|_{\rm F}^{2}$ aims to
control and $s^{\prime}(\lambda)$ measures the sensitivity of the target
quantity $s(\lambda)$ with respect to $\lambda$. Thus, our tuning method
chooses $\lambda$ corresponding to a prespecified sensitivity level $\delta$.

To illustrate the advantage of tuning $\lambda$ via $s^{\prime-1}(\cdot)$, we
implement RQMF with 50 different $\lambda$ evenly spaced in $[0,0.1]$. We use
samples drawn from the sine curve, that is, $(t_{i},\sin(t_{i}))+\epsilon_{i}$
with $\\{t_{i}\\}_{i=1}^{21}$ evenly distributed in
$[\frac{\pi}{3},\frac{2\pi}{3}]$ and $\epsilon_{i}\overset{\rm
i.i.d.}{\sim}{\cal N}({\bf 0},0.03^{2}I)$. For each $\lambda$, we implement
RQMF and compute the error, that is, the average distance between the fitted
data points and the underlying truth. Also, we calculate the values of
$s(\lambda)$ and $\delta(\lambda)=-s^{\prime}(\lambda)$. Figure 4 displays the
error against different $\lambda$, $s(\lambda)$, and $\delta(\lambda)$. It
shows that the error versus $\delta(\lambda)$ curve is the flattest near the
optimal error. To achieve a prespecified error, say $1.2\times 10^{-3}$, the
feasible choice of $\delta(\lambda)$ has a much wider range than that of
$\lambda$ or $s(\lambda)$. Thus, it is easier to achieve good performances of
RQMF by choosing $\delta(\lambda)$ rather than $\lambda$ or $s(\lambda)$. In
practice, we choose $\delta>0$ as a constant smaller than $-s^{\prime}(0)$,
where $\Phi$ determining the function $s(\cdot)$ is given by LMF.

######  Proposition 6.

Suppose $\Phi$ is fixed with $\Phi\Phi^{T}=I_{d}$ and $\Phi{\bf 1}_{m}=0$.
Define $\widetilde{R}(\lambda)$ by (32) and
$s(\lambda)=\|\widetilde{R}(\lambda)J\|_{\rm F}^{2}$. Then
$s^{\prime}(\lambda)\leq 0$ and $s^{\prime\prime}(\lambda)\geq 0$.
Furthermore, the strict inequalities $s^{\prime}(\lambda)<0$ and
$s^{\prime\prime}(\lambda)>0$ hold if $s(\lambda)>0$ and $\lambda>0$.

The proof of Proposition 6 is collected in the Appendix.

##  5 Applications to Manifold Learning

This section applies the RQMF algorithm to manifold learning problems. Assume
data $\\{x_{i}\\}_{i=1}^{m}\subseteq\mathbb{R}^{D}$ are generated near an
unknown smooth manifold $\cal M$ of intrinsic dimension $d$. Here we no longer
assume all data belong to the same local chart. Instead, we assume data belong
to a union of several local charts. To recover the underlying manifold, we can
apply the RQMF algorithm for each local chart.

The performance of this divide-and-conquer strategy depends on choices of
specific local charts. Besides, the quality of the fitted points on a single
chart cannot be guaranteed uniformly: recovering the central region tends to
be of higher quality than recovering the marginal region. Also, for a data
point belonging to multiple local charts, the fitted points in different
charts are different and it is hard to determine which one is the best. To
address these challenges, we propose an improved divide-and-conquer strategy.
This strategy constructs a local chart for each data point $y$ by finding its
nearest $K$ samples or by $\mathcal{N}(y,a)=\\{x_{i}\mid\|x_{i}-y\|\leq a\\}$
for some $a>0$. Then for each target sample, we denoise this particular data
point by applying the RQMF algorithm to the corresponding chart. Compared with
the original divide-and-conquer strategy, our strategy treats each data point
as an individual problem and improves the accuracy.

In a more general form, we could use a kernel function $K_{h}(\cdot,\cdot)$ to
assign a closer point with higher importance, where $h$ is the bandwidth. For
a target sample $y$, we modify the loss function in (30) as

$\displaystyle\min_{R,\Phi\atop\Phi\Phi^{T}=I,\Phi{\bf 1}_{m}^{T}={\bf
0}}\ell_{\lambda,y,h}(R,\Phi)=$
$\displaystyle\|(X-RT(\Phi))W_{h}^{1/2}(y)\|_{\rm F}^{2}$ (34)
$\displaystyle+\lambda\|RJ\|_{\rm F}^{2},$

where $X=(x_{1},\ldots,x_{m})$ is the global data matrix and
$W_{h}^{1/2}(y)\in\mathbb{R}^{m\times m}$ is a diagonal weight matrix with the
$i$-th diagonal element equal to $K_{h}^{1/2}(x_{i},y)$. If we choose
$K_{h}(x,y)=1_{\|x-y\|\leq a}$, then (34) reduces to the improved divide-and-
conquer strategy mentioned above. It is also possible to use other kernel
functions, such as the Gaussian kernel. To distinguish the kernel RQMF model
from the previous equal-weight RQMF, we use RQMF-E to represent equal-weigth
RQMF and RQMF-K to represent RQMF with weights determined by a kernel.

We also discuss how to tune $\lambda$ for each sub-problem (34). Picking the
same $\lambda$ for all sub-problems is not desirable since the best $\lambda$
depends on the weights $W_{h}^{1/2}(y)$ and the curvature of the underlying
truth, which vary as $y$ changes. Instead, we suggest using the tuning method
proposed in Section 4.1, which picks $\lambda=s^{\prime-1}(-\delta)$ for the
same prescribed sensitivity level $\delta>0$ for all sub-problems. The same
$\delta$ would result in different $\lambda$’s for different charts, and this
strategy often leads to better fitting accuracies in our experience.

TABLE I: Comparisons of different methods on the synthetic spherical dataset
in terms of MSE and SD (in bracket) with varying $K$.

$K$ | 7 | 10 | 13 | 16 | 19 | 22 | 25 | 28
---|---|---|---|---|---|---|---|---
RQMF-E |  0.0243 (0.0325) | 0.0165 (0.0227) |  0.0122 (0.0172) |  0.0115 (0.0164) |  0.0148 (0.0256) |  0.0130 (0.0222) |  0.0149 (0.0344) |  0.0156 (0.0356)
RQMF-K |  0.0188 (0.0281) |  0.0170 (0.0270) |  0.0155 (0.0231) |  0.0148 (0.0233) |  0.0153 (0.0242) |  0.0159 (0.0252) | 0.0170 (0.0266) |  0.0185 (0.0280)
Local PCA | 0.0437 (0.0521) | 0.0437 (0.0522) | 0.0435 (0.0520) | 0.0432 (0.0519) | 0.0434 (0.0521) | 0.0434 (0.0518) | 0.0434 (0.0517) | 0.0434 (0.0517)
KDE | 0.0333 (0.0480) | 0.0298 (0.0469) | 0.0302 (0.0483) | 0.0307 (0.0482) | 0.0323 (0.0493) | 0.0342 (0.0499) | 0.0369 (0.0492) | 0.0389 (0.0501)
LOG-KDE | 0.0278 (0.0417) | 0.0192 (0.0350) | 0.0159 (0.0344) | 0.0155 (0.0348) | 0.0168 (0.0349) | 0.0192 (0.0370) | 0.0230 (0.0395) | 0.0275 (0.0415)
Mfit | 0.0392 (0.0463) | 0.0333 (0.0424) | 0.0262 (0.0347) | 0.0215 (0.0314) | 0.0183 (0.0302) | 0.0160 (0.0299) |  0.0154 (0.0303) |  0.0185 (0.0376)
Moving LS | 0.0420 (0.0844) | 0.0673 (0.1155) | 0.1017 (0.1474) | 0.1506 (0.1666) | 0.1898 (0.1819) | 0.2264 (0.2033) | 0.2544 (0.2060) | 0.2642 (0.2072)
Numbers in bold and underlined are the best and second-best results for each
column’s setting, respectively.

##  6 Numerical Experiments

This section presents numerical experiments on a synthetic manifold learning
dataset and two real datasets, including the MNIST handwritten dataset and a
cryogenic electron microscopy dataset, to examine the finite-sample
performance of the proposed method. Our goal is to reconstruct the underlying
manifold from noisy data and compare our method with five commonly used
methods. We first briefly describe these five competitors.

  * •

Local PCA For any $x$, it first finds the $K$ nearest data points
$\\{x_{i_{1}},...,x_{i_{K}}\\}$ and then computes the covariance matrix
$M=\frac{1}{K}\sum_{k=1}^{K}(x_{i_{k}}-c_{x})(x_{i_{k}}-c_{x})^{T}\in\mathbb{R}^{D\times
D}$, where $c_{x}=\sum_{k=1}^{K}x_{i_{k}}/K$ is the center of these samples.
Denote by $P\in\mathbb{R}^{D\times D}$ the projection matrix corresponding to
the space spanned by the $d$ principle eigenvectors of $M$. The denoised point
of $x$, a point on the estimated manifold “projected” from $x$, is given by
$x_{\textnormal{new}}=c_{x}+P(x-c_{x})$.

  * •

KDE & LOG-KDE Ridge Estimation Both methods are special cases of the
nonparametric ridge estimation method [21]. Let
$\hat{p}(x)=\sum_{i}K_{h}(x,x_{i})$ be the kernel density estimation (KDE)
with the kernel function $K_{h}(\cdot,\cdot)$ and the bandwidth $h$. KDE ridge
estimator estimates the ridge:

$\displaystyle{\rm ridge}\coloneqq\\{x\mid$
$\displaystyle\Pi^{\perp}(\nabla^{2}\hat{p}(x))\nabla\hat{p}(x)={\bf 0},$ (35)
$\displaystyle\lambda_{d+1}(\nabla^{2}\hat{p}(x))<0\\},$

where $\Pi^{\perp}(\nabla^{2}\hat{p}(x))=I-UU^{T}$ with
$U\in\mathbb{R}^{D\times d}$ given by the top $p$ principal eigenvectors of
$\nabla^{2}\hat{p}(x)$. Although the ridge in (35) does not admit close-form
solutions, we may use the subspace constrained mean shift (SCMS) algorithm to
find the denoised point of any point $x$ and thus the ridge [20]. Similarly,
the LOG-KDE ridge estimation is merely KDE ridge estimation with $\hat{p}(x)$
replaced by $\log\hat{p}(x)$.

  * •

Mfit Mfit, proposed by [22], estimates the following manifold:

$\displaystyle\left\\{x\mid\Pi_{x}\big{(}\sum_{i}\alpha(x,x_{i})\Pi_{i}(x-x_{i})\big{)}={\bf
0}\right\\},$

where $\alpha(\cdot,\cdot)$ is a weight function, $\Pi_{i}$ is the projection
matrix onto the approximate normal space at $x_{i}$, and $\Pi_{x}$ is the
projection matrix corresponding to the top $D-d$ principal eigenvectors of the
matrix $\sum_{i}\alpha(x,x_{i})\Pi_{i}$. Again, we may use the SCMS algorithm
to solve this problem.

  * •

Moving LS The moving least square (LS) consists of two steps [27]. For any
$x\in\mathbb{R}^{D}$, we first find the $d$-dimensional hyperplane $\cal H$ in
$\mathbb{R}^{D}$ minimizing the following quantity

$\displaystyle{\cal L}_{1}({\cal H})=\min_{q\in{\cal H},x-q\perp{\cal
H}}\sum_{i}\alpha(q,x_{i})\rho^{2}(x_{i},{\cal H}),$

where $\alpha(\cdot,\cdot)$ is a weight function and $\rho(x_{i},{\cal H})$ is
the distance between $x_{i}$ and ${\cal H}$. We can construct a coordinate
system on $\cal H$ with origin $q$, where $q$ is the projection point of $x$
on $\cal H$. Using this coordinate system, we can obtain the $d$-dimensional
configuration $x_{i}^{\prime}$ of $x_{i}$ by projecting $x_{i}$ to $\cal H$.
Second, we fit a polynomial function $p:\mathbb{R}^{d}\to\mathbb{R}^{D}$ of a
given degree by minimizing the following weighted squares

$\displaystyle{\cal
L}_{2}(p)=\sum_{i}\alpha(q,x_{i})\|p(x_{i}^{\prime})-x_{i}\|_{2}^{2}.$

The denoised point of $x$ is then given by $p(0)$. In our experiments, we fix
the degree of $p$ as two.

Figure 5: An illustration of the impact of $\delta$ and $K$ for spherical
fitting for RQMF-E (left) and RQMF-K(right).  Figure 6: An illustration of
local fitted surface under the impact of $\delta$ for RQMF-E when $K=18$.

###  6.1 A Synthetic Example

In this subsection, we compare RQMF-E and RQMF-K with the above five
competitors in a synthetic spherical fitting experiment. We simulate the noisy
data $\\{x_{i}\\}_{i=1}^{240}$ by generating 240 points uniformly from the
unit sphere $\cal S$ in $\mathbb{R}^{3}$ first and then adding independent
noises following ${\cal N}(0,\sigma^{2}I)$ with $\sigma=0.2$. All algorithms
take in the noisy data $\\{x_{i}\\}$ and then output the denoised data
$\\{\widehat{x}_{i}\\}$. To measure the performance of different algorithms,
we use the mean squared error (MSE) and standard derivation (SD):

$\displaystyle{\rm MSE}=\sum_{i=1}^{m}\|\widehat{x}_{i}-P_{\cal
S}(\widehat{x}_{i})\|^{2}_{2}/m,$ $\displaystyle{\rm
SD}=\sqrt{\frac{1}{m}\sum_{i=1}^{m}(\|\widehat{x}_{i}-P_{\cal
S}(\widehat{x}_{i})\|_{2}^{2}-{\rm MSE})^{2}},$

where $P_{\cal S}(\cdot)$ is the projector onto the sphere.

Figure 7: The first row displays 12 examples from the original MNIST dataset. The second to the sixth rows collect the results for the RQMF-K, KDE, LOG-KDE, Mfit, and Moving LS algorithms, respectively. TABLE II: Comparison of the smoothness measured by the average of the nonzero of the absolute value of $I*w$ corresponding to 12 example images  Image ID | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---|---
Original | 1.0637 | 1.2264 | 1.0306 | 0.8760 | 1.0332 | 0.9499 | 1.0246 | 1.0261 | 1.0068 | 1.0267 | 0.7178 | 1.0326
RQMF-K | 0.4848 | 0.5795 | 0.5168 | 0.5785 | 0.5440 | 0.4851 | 0.5733 | 0.5319 | 0.6011 | 0.5388 | 0.5327 | 0.5470
KDE | 1.0875 | 0.8876 | 0.7574 | 0.9382 | 0.8387 | 0.9246 | 0.7489 | 0.6942 | 0.6687 | 0.9070 | 0.6320 | 0.7359
LOG-KDE | 1.1395 | 0.8861 | 0.7732 | 0.9460 | 0.8606 | 0.9329 | 0.7650 | 0.7523 | 0.6979 | 0.8812 | 0.6361 | 0.7764
Mfit | 1.1769 | 0.9158 | 0.8567 | 0.9026 | 0.9520 | 0.9705 | 0.7956 | 0.8274 | 0.9107 | 0.8562 | 0.6189 | 0.7976
Moving LS | 1.1926 | 1.3185 | 1.4126 | 1.0072 | 1.2117 | 1.1529 | 1.4050 | 1.3821 | 1.4752 | 1.2451 | 0.9490 | 1.4394
Numbers in bold and underlined are the best and second-best results for each
column’s setting, respectively.

We briefly discuss how RQMF-E and RQMF-K are implemented. RQMF-E denoises each
data point $y$ using the $K$ nearest neighbors of $y$, where $K$ is a tuning
parameter. To avoid overfitting, we require $K$ to be larger than the rank
$(d^{2}+3d+2)/2$ of $RT(\Phi)$. For RQMF-K, we use the Gaussian kernel and set
the bandwidth as $h=d_{K}/3+3$, where $d_{K}$ is the distance from $y$ to its
$K$-th nearest neighbor. For both RQMF algorithms, we set the regularization
parameter $\lambda$ as $\lambda=s^{\prime-1}(-\delta)$, where $\delta$ is a
tuning parameter. To select $K$ and $\delta$ for both RQMFs, we compute the
MSEs of both RQMF-E and RQMF-K with different $(K,\delta)$ as shown in Figure
5. When $K$ is small, the local data approximates a plane and thus a smaller
$\delta$ (larger $\lambda$) yields a better performance for RQMF-E. When $K$
is large, the local data exhibits nonlinear structures, thus it is better to
use a larger $\delta$ (smaller $\lambda$) for RQMF-E. The effect of $\delta$
is visualized in Figure 6. It shows that smaller $\delta$ tends to fit flatter
planes in comparison with larger $\delta$, which coincides with the phenomenon
in Figure 5. To characterize the good performance region for $\delta$ and $K$
in Figure 5, we choose $\delta=\max\\{1,8K-125\\}$ to determine $\delta$ based
on $K$ for RQMF-E in this experiment. On the other hand, the performance of
RQMF-K is relatively robust to the choice of $\delta$, thus we fix
$\delta=100$ for different choices of $K$.

Now we compare RQMF-E and RQMF-K with their competitors. The results are
collected in Table I. The results indicate that RQMF-E and RQMF-K outperform
other methods for a wide range of $K$. When $K=7$, RQMF-K achieves the best
performance and when $10<K<28$, RQMF-E achieves the best performance among all
methods. If we focus on the best performance of different algorithms, the
RQMF-E is still favored with the minimal ${\rm MSE}=0.0115$ when $K=16$. The
superior performances of RQMF demonstrate the benefits of using the curvature
information in the denosing procedure. It is also worth noting that RQMF
outperforms Moving LS, which also fits a quadratic polynomial in its second
step. This is possibly due to fact that RQMF iteratively updates the local
representations $\Phi$ of data points, while Moving LS only uses the local
coordinates learned in its first step. The estimation error of the local
coordinates learned in Moving LS could lead to a degradation of the final
fitting accuracy.

###  6.2 An Application to the MNIST Handwritten Digit Dataset

This subsection compares RQMF-K and its competitors on the MNIST handwritten
digit dataset [28]. Each image in the dataset consists of $28\times 28$
pixels. We use $g(a,b)$ to denote the grey value of an image at pixel $(a,b)$
and each image is determined by such a function $g$. Only pixels with nonzero
grey values are considered, so the dataset for each image is given by
$\\{x_{i}=(a_{i},b_{i})\in\mathbb{R}^{2}\mid g(a_{i},b_{i})>0\\}$. In this
way, each image can be viewed as a perturbed one-dimensional manifold in
$\mathbb{R}^{2}$ and our goal is to recover the underlying manifold, which is
also referred to as the principal curve [20].

The pixel closer to the principal curve tends to have a larger grey value.
Thus, it is natural to use the grey values as weights in (34). Specifically,
for each $y$, we set the diagonal weight matrix $W_{h}(y)$ in (34) by
$(W_{h}(y))_{ii}=K_{h}(y,x_{i})g(x_{i})/s$ for some constant $s$, where
$K_{h}(\cdot,\cdot)$ is a Gaussian kernel and the bandwidth $h$ is given by
the distance of $y$ to $y$’s $K$-th nearest pixel. We use $\delta=100$ to tune
the parameter $\lambda=s^{\prime-1}(\delta)$.

Since the smoooth curve contains the most significant signal and the isolated
points can be thought as noise, we measure the smoothness of the image using
the convolution of the original image with the Laplace operator

$\displaystyle w=\left[\begin{array}[]{ccc}0&1&0\\\ 1&-4&1\\\
0&1&0\end{array}\right].$

The obtained matrix $I*w$ is a discrete version of the Laplace operator
defined for function $f$, i.e.,
$\Delta^{2}f(x,y)=\frac{\partial^{2}f}{\partial
x^{2}}+\frac{\partial^{2}f}{\partial y^{2}}$. The average value of $I*w$
represents the degree of smoothness of the image $I$. We report the average of
the nonzero of the absolute value in $I*w$ in Table II.

Figure 8: Visualization the fitting and denoised results of 16 randomly chosen
images by RQMF-E and five related manifold learning methods with $K=60$.

For each image, we apply all six manifold learning algorithms to recover the
principal curve, and the results are displayed in Figure 7 and Table II. It
turns out that the denoised images by RQMF-K are smoother than the original
images and images output by KDE, LOG-KDE, Mfit. While Moving LS also produces
smoother images, it tends to twist the original images too much; see images of
6,8,9. In contrast, RQMF preserves the major contour of the original images
much better.

###  6.3 An Application to Cryo-EM

This subsection compares RQMF-E and RQMF-K with its competitors on the Cryo-EM
dataset [29]. This dataset consists of $n=2000$ images with shape $64\times
64$. Each image is modeled as a vector in $\mathbb{R}^{4096}$ with elements
given by the grey values on all 4096 pixels. The whole dataset is then
represented as $\\{\iota_{i}\\}_{i=1}^{m}\subseteq\mathbb{R}^{4096}$, where
$\iota_{i}$ denotes the $i$-th image. The dataset inherently resides on a
lower-dimensional manifold in $\mathbb{R}^{4096}$ due to the image generation
and processing of Cryo-EM, such as rotation, projection, and blurring by
convolution [22]. In our experiment, we take the original data as the
underlying manifolds, add noises to the original data, apply all six manifold
learning algorithms to the noisy data, and finally compare their recovery
accuracies.

Figure 9: The performance of RQMF-E with different $K$ and $\delta$ for Cryo-EM.  TABLE III:  Comparisons of six manifold learning algorithms on the Cryo-EM dataset. We report MSE and SD (in bracket) with varying $K$. $K$ | 40 | 44 | 48 | 52 | 56 | 60
---|---|---|---|---|---|---
RQMF-E |  1.5757(0.7065) |  1.4203(0.6669) |  1.3482(0.6013) |  1.3791(0.6363) |  1.4134(0.6465) |  1.3764(0.6623)
RQMF-K | 2.1407(1.1057) | 2.1008(1.0804) | 2.0718(1.0612) | 2.0310(1.0385) | 1.9889(1.0064) | 1.9399(0.9618)
Local-PCA | 1.8296(0.7531) | 1.8091(0.7211) | 1.7364(0.7662) | 1.7423(0.7476) | 1.7822(0.7958) | 1.8130(0.7553)
KDE | 4.5179(1.5050) | 4.4408(1.4786) | 4.3820(1.4143) | 4.3258(1.4975) | 4.3213(1.4728) | 4.2659(1.4698)
LOG-KDE | 1.9323(0.9683) | 1.8667(0.9296) | 1.8461(0.9188) | 1.8018(0.8514) | 1.7974(0.8288) | 1.7955(0.8242)
Mfit | 2.4032(1.1298) | 2.3861(1.1216) | 2.3477(1.0988) | 2.2902(1.1229) | 2.2852(1.0957) | 2.2225(1.0391)
Moving-LS |  1.5518(1.0340) |  1.4685(0.7865) |  1.5185(0.9541) |  1.5515(1.2535) |  1.5083(1.1097) |  1.4567(0.9147)
Numbers in bold and underlined are the best and second-best results for each
column’s setting, respectively.

It is computationally expensive to directly fitting the manifold in
$\mathbb{R}^{4096}$ using any manifold learning algorithm. To reduce the
dimensionality, we approximate the original dataset by a $D$-dimensional
subspace such that $\iota_{i}\approx Ux_{i}$, where
$U\in\mathbb{R}^{4096\times D}$ denotes $D$ principal eigenvectors of
$S=\frac{1}{m}\sum_{i}\iota_{i}\iota_{i}^{T}$ and
$x_{i}=U^{\top}\iota_{i}\in\mathbb{R}^{D}$. We could fit a $d$-dimensional
manifold in $\mathbb{R}^{D}$ and map such a manifold to the pixel space
$\mathbb{R}^{4096}$ via $U$. In what follows, we fix $D=20$. We construct the
noisy dataset as $\\{\iota^{\prime}_{i}=Ux^{\prime}_{i}\\}_{i=1}^{m}$ with
$x^{\prime}_{i}$ given by

$\displaystyle x_{i}^{\prime}=x_{i}+\epsilon_{i},\quad\epsilon_{i}\sim
N(0,\sigma^{2}I_{D}),\quad i=1,\ldots,m.$

Next, we recover a $5$-dimensional manifold in $\mathbb{R}^{D}$ by applying
all six manifold learning algorithms to $\\{x^{\prime}_{i}\\}_{i=1}^{m}$. For
each method and each sample, we use the nearest $K$ samples to fit the local
manifold. In addition, we use the Gaussian kernel in (34) and for each $y$, we
set the bandwidth $h$ as $\|y-x_{i_{K}}\|_{2}$, where $x_{i_{K}}$ is the
$K$-th nearest neighbour of $y$. Figure 8 visualizes 16 denoised images using
these manifold learning methods with $K=60$.

To evaluate the performance, we use the following mean squared error and
standard deviation of the error between the solution and the real image:

$\displaystyle{\rm
MSE}=\frac{1}{m}\sum_{i=1}^{m}\|\iota_{i}-\widehat{\iota}_{i}\|_{2}^{2},$
$\displaystyle{\rm
SD}=\sqrt{\frac{1}{m}\sum_{i=1}^{m}(\|\iota_{i}-\widehat{\iota}_{i}\|_{2}^{2}-{\rm
MSE})^{2}},$

where $\hat{\iota}_{i}=U\hat{x}_{i}$ for all $i$ and $\hat{x}_{i}$ is the
$i$-th fitted point in $\mathbb{R}^{D}$. Figure 9 displays the MSE of RQMF-E
with different $K$ and $\delta$. It shows that RQMF-E achieves the best
performance when $K=48$ and $\delta=50$.

To compare different methods, we collect the MSEs of all methods with varying
$K$ in Table III. It can be seen that RQMF-E outperforms other methods in
terms of MSE and SD in most settings. If we focus on the best performance of
each method, RQMF-E is again favored with an error 1.3482 when $K=48$.
Therefore, by taking the curvature information into account, RQMF-E exhibits a
stronger expressive ability than other methods and thus achieves better
denoising performance on this dataset.

##  7 Concluding Remarks

This paper proposes a quadratic matrix factorization framework to learn the
structure of the observed data. We develop an alternating minimization
algorithm to solve the non-convex quadratic matrix factorization problem as
well as a regularized version. Theoretical convergence properties are
established. We also present a novel transformation-based parameter-tuning
method for regularized quadratic matrix factorization and intuitively argue
its advantages over naively tuning the original regularization parameter.
Furthermore, we apply the proposed methods to manifold learning problems. We
demonstrate the superiority of the proposed method numerically in a synthetic
manifold learning dataset and two real datasets, i.e., the MNIST handwritten
dataset and a cryogenic electron microscopy dataset.

There are several interesting directions for future research. First, our work
and most related works assume the intrinsic dimension $d$ is known a priori,
while this information is often not available in practice. Thus, it remains an
important question to estimate the intrinsic dimensionality $d$ under the
quadratic matrix factorization framework. It would also be interesting to
characterize the impact if $d$ is misspecified. Second, the noises in the
signal-plus-noise model may be heavy-tailed or even adversarial, so it is
important to develop robust algorithms. Third, non-negative constraints are
widely used in linear matrix factorization [6, 12]. It is interesting to study
how non-negative constraints can be used in QMF to enhance performance.

## References

  * [1] Pentti Paatero and Unto Tapper.  Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values.  Environmetrics, 5(2):111–126, 1994.
  * [2] Jialu Liu, Chi Wang, Jing Gao, and Jiawei Han.  Multi-view clustering via joint nonnegative matrix factorization.  In Proceedings of the 2013 SIAM international conference on data mining, pages 252–260. SIAM, 2013.
  * [3] Qi Wang, Xiang He, Xu Jiang, and Xuelong Li.  Robust bi-stochastic graph regularized matrix factorization for data clustering.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(1):390–403, 2022.
  * [4] Yehuda Koren, Robert Bell, and Chris Volinsky.  Matrix factorization techniques for recommender systems.  Computer, 42(8):30–37, 2009.
  * [5] William L Hamilton, Rex Ying, and Jure Leskovec.  Representation learning on graphs: Methods and applications.  arXiv preprint arXiv:1709.05584, 2017.
  * [6] Daniel Lee and H Sebastian Seung.  Algorithms for non-negative matrix factorization.  Advances in neural information processing systems, 13, 2000.
  * [7] Chris HQ Ding, Tao Li, and Michael I Jordan.  Convex and semi-nonnegative matrix factorizations.  IEEE transactions on pattern analysis and machine intelligence, 32(1):45–55, 2008.
  * [8] Chris Ding, Xiaofeng He, and Horst D Simon.  On the equivalence of nonnegative matrix factorization and spectral clustering.  In Proceedings of the 2005 SIAM international conference on data mining, pages 606–610. SIAM, 2005.
  * [9] Eddie Aamari and Clément Levrard.  Nonasymptotic rates for manifold, tangent space and curvature estimation.  The Annals of Statistics, 47(1):177–204, 2019.
  * [10] Yariv Aizenbud and Barak Sober.  Non-parametric estimation of manifolds from noisy data.  arXiv preprint arXiv:2105.04754, 2021.
  * [11] Yu-Xiong Wang and Yu-Jin Zhang.  Nonnegative matrix factorization: A comprehensive review.  IEEE Transactions on knowledge and data engineering, 25(6):1336–1353, 2012.
  * [12] Nicolas Gillis and Abhishek Kumar.  Exact and heuristic algorithms for semi-nonnegative matrix factorization.  SIAM Journal on Matrix Analysis and Applications, 36(4):1404–1424, 2015.
  * [13] Maryam Abdolali and Nicolas Gillis.  Simplex-structured matrix factorization: Sparsity-based identifiability and provably correct algorithms.  SIAM Journal on Mathematics of Data Science, 3(2):593–623, 2021\.
  * [14] Zhirong Yang and Erkki Oja.  Unified development of multiplicative algorithms for linear and quadratic nonnegative matrix factorization.  IEEE Transactions on Neural Networks, 22(12):1878–1891, 2011.
  * [15] Zhirong Yang and Erkki Oja.  Quadratic nonnegative matrix factorization.  Pattern Recognition, 45(4):1500–1510, 2012.
  * [16] Sam T Roweis and Lawrence K Saul.  Nonlinear dimensionality reduction by locally linear embedding.  science, 290(5500):2323–2326, 2000.
  * [17] Joshua B Tenenbaum, Vin De Silva, and John C Langford.  A global geometric framework for nonlinear dimensionality reduction.  science, 290(5500):2319–2323, 2000.
  * [18] Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  Neural computation, 15(6):1373–1396, 2003.
  * [19] Ronald R Coifman and Stéphane Lafon.  Diffusion maps.  Applied and computational harmonic analysis, 21(1):5–30, 2006.
  * [20] Umut Ozertem and Deniz Erdogmus.  Locally defined principal curves and surfaces.  Journal of Machine learning research, 12(Apr):1249–1286, 2011.
  * [21] Christopher R Genovese, Marco Perone-Pacifico, Isabella Verdinelli, Larry Wasserman, et al.  Nonparametric ridge estimation.  The Annals of Statistics, 42(4):1511–1545, 2014.
  * [22] Charles Fefferman, Sergei Ivanov, Yaroslav Kurylev, Matti Lassas, and Hariharan Narayanan.  Fitting a putative manifold to noisy data.  In Conference On Learning Theory, pages 688–720, 2018.
  * [23] Dian Gong, Fei Sha, and Gérard Medioni.  Locally linear denoising on image manifolds.  In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 265–272. JMLR Workshop and Conference Proceedings, 2010.
  * [24] Yen-Chi Chen, Christopher R Genovese, Larry Wasserman, et al.  Asymptotic theory for density ridges.  The Annals of Statistics, 43(5):1896–1928, 2015.
  * [25] Alan Edelman, Tomás A Arias, and Steven T Smith.  The geometry of algorithms with orthogonality constraints.  SIAM journal on Matrix Analysis and Applications, 20(2):303–353, 1998.
  * [26] Amir Beck.  First-order methods in optimization.  SIAM, 2017.
  * [27] Barak Sober and David Levin.  Manifold approximation by moving least-squares projection (mmls).  Constructive Approximation, 52(3):433–478, 2020.
  * [28] Li Deng.  The mnist database of handwritten digit images for machine learning research.  IEEE signal processing magazine, 29(6):141–142, 2012.
  * [29] Xiao-Chen Bai, Greg McMullan, and Sjors HW Scheres.  How cryo-em is revolutionizing structural biology.  Trends in biochemical sciences, 40(1):49–57, 2015.

Zheng Zhai received his Ph.D. degree from the School of Mathematical Sciences,
Zhejiang University, and he is currently working as a postdoc researcher at
the Department of Statistical Sciences, University of Toronto. His current
research interests include unsupervised learning and manifold learning.
---
Hengchao Chen is currently pursuing his Ph.D. degree in statistics at the
University of Toronto. His research interests include matrix factorization and
manifold learning.
---
Qiang Sun received his Ph.D. from the University of North Carolina at Chapel
Hill. He was an associate research scholar before he joined the University of
Toronto as a faculty member. He is currently an associate professor of
statistics at the University of Toronto. His research lies in the intersection
of statistics, optimization, and learning.
---


# Comment on: In vitro prediction of the lower/upper-critical biofluid flow
choking index and in vivo demonstration of flow choking in the stenosis artery
of the animal with air embolism

T. Podgorski<EMAIL_ADDRESS>Université Grenoble
Alpes, CNRS, Grenoble INP, LRP, F-38000 Grenoble, France

###### Abstract

Sanal Kumar et al. Sanal Kumar  _et al._ (2022a) present a model of transonic
compressible flows based on ideal gas theory that is irrelevant to biofluid
flow and there are flaws in the general reasoning. In addition, the
experimental attempts do not show any evidence of supersonic flow and do not
provide any support for the flawed theory.

††preprint: AIP/123-QED

In a recent paper, Sanal Kumar et al. Sanal Kumar  _et al._ (2022a) propose
predictions of a flow choking index that is supposed to be relevant to the
understanding of various vascular strokes and diseases which in their opinion
are supposed to be triggered when a critical Mach number is reached. They also
provide their own interpretation of issues related to air embolisms and
describe experimental and numerical results which are claimed to be supporting
their theory. However, besides the fact that the logical structure of the
article is fuzzy, there are several major issues with this work, from the
framework and applicability of the model, to logical fallacies in its
development and conclusions as well as false statements on blood properties
and irrelevant or non-physical experimental results. Only the main and most
salient controversial points are discussed here. Note that some of these
criticisms also hold for their previous paper in this journalSanal Kumar  _et
al._ (2022b) as it is based on the exact same erroneous model.

##  I Relevance of the model

The main model equations (Eqs (1-5)) which were derived in one of the authors’
previous publications Kumar  _et al._ (2018) are simply based on a mass
conservation equation and the classic equation of state and expression of the
speed of sound for an ideal gas. This assumption on which the model is based
is already a major flaw since liquids, biofluids including blood can certainly
not be considered as ideal gases (if it were so, following the well known
$PV=nRT$ law, our body volume would double when going from sea level to
Himalayan mountains).

In other words, even considering that (in the strict sense) any material,
solid or fluid, is compressible, the ideal gas equation is not even remotely
an approximation of their behavior. Blood being mainly composed of water, its
compressibility is close to that of water, a well-tabulated property which is
around $4.4\times 10^{-10}$ Pa-1 at body temperature and around atmospheric
pressureFine and Millero (1973), with only little variations with pressure and
temperature. This is $10^{5}$ times lower that the compressibility of ideal
gases at atmospheric pressure. This fully justifies that blood is safely
considered nearly incompressible for all practical flow conditions with solid
proof from a plethora of theoretical, numerical and experimental studies
available in the literature. Equations (1-5) are therefore irrelevant to
biofluids and only apply to ideal gases.

Now, even supposing that these equations could be modified by replacing the
ideal gas law by the equation of state for a real fluid, the general approach
which consists in looking for critical Mach numbers (i.e. suspecting that
supersonic flows occur in the cardiovascular system) is bound to fail because
the well documented range of blood flow velocities is between 1 mm/s
(capillaries) and less than 1 m/s (arteries) everywhere in the circulatory
systemRobertson, Sequeira, and Kameneva (2008); Podgorski (2022) while on the
other hand the speed of sound in blood, which is also well documented in
textbooksDuck (1990) and widely used in medical ultrasonic techniques is close
to that of water, that is about 1500 m/s. In other words, the Mach number in
blood flows is always less than $10^{-3}$, meaning that the "Sanal flow
choking" theory presented in the commented paper (and other similar papers by
the same group of authors), even if it were to be corrected for real fluids,
is irrelevant to cardiovascular flows.

Following this clarification, it is clear that none of the in silico results
presented in Figures 7-11 of the paper represent realistic blood flow in
arteries. No numerical simulation of blood flow in arteries using relevant
values of physical parameters (viscosity, density) and boundary conditions
(flow rates, pressures, velocities, channel dimensions), even including
turbulence modeling Lui  _et al._ (2020) would lead to the flow patterns shown
in Figures 7-11 where the Mach number goes up to 2. Unfortunately, none of the
parameter values, boundary conditions and length scales used in these
simulations are given in the article (the axial distance axis has no units in
any of these figures).Only the video linked to Fig. 8 shows a scale bar, from
which we note that their simulated artery has an irrelevant diameter between
10 cm and 20 cm. In addition, the inlet Mach number in these figures also has
completely unrealistic values of order 1, meaning that the authors assume that
blood velocity has nonsensical values of the order of kilometers per second
when entering the depicted artery section. Stunningly, some of these figures
are the exact same as the ones they show in other papers about aircraft or
rocket aerodynamics, for instance Figure 7, which according to the caption
shows numerical simulations of blood flow in an artery, is exactly the same as
Figure 1 in a previous paper on aircraft dynamicsSanal Kumar  _et al._ (2021),
and has also been used in another previous paper where it is supposed to apply
to different casesSanal Kumar  _et al._ (2022b). Therefore, these simulation
results which were performed with unrealistic boundary conditions and
unrealistic parameter values must be totally discarded with regard to their
applicability to blood flows.

As an additional note, model equations are based on circular reasoningKumar
_et al._ (2018): they assume that a boundary layer develops in the channel
(through unspecified mechanisms), which effectively reduces the cross-section
of the flow, which because of mass conservation leads to an increase of the
flow velocity, leading to equation 5 which is nothing more than a
reformulation of the trivial relationship between the local hydraulic diameter
and local velocity in a channel imposed by mass conservation. This triviality
can obviously not predict the existence nor the thickness of this boundary
layer as a function of distance along the flow axis since it is a single
equation relating two unknowns (the so-called boundary layer thickness and the
local average velocity). The specific heat capacities of water (blood’s main
constituent) are $C_{p}=4.18$ kJ/(kg.K) and $C_{p}=4.09$ kJ/(kg.K) at
$T=37^{\circ}$C as can be found in any handbookRažnjević (1995) (values are
similar for bloodMcintosh and Anderson (2010)), the heat capacity ratio is
then $\gamma=1.02$. It is indeed always close to 1 in all physiologically
relevant cases contrary to what the authors want us to believe in the
introduction (it only goes up to $1.03$ at $T=42^{\circ}$C, lethal
temperature). Therefore the trivial mass conservation equation (Eq. 5a) almost
reduces to $Ud^{2}\simeq$ constant and no "choking" or supersonic velocity is
expected with realistic inlet velocities, pressures and channel diameters,
even in the most pathological cases. Finally, equation 1 cannot represent
systolic to diastolic pressure ratios since their modeling does not involve
any time dependency or transient dynamics but only steady flow: it is an
illicit use and combination of classic equations of compressible fluid theory
taken out of context.

##  II Flawed experiments

In an attempt to support their theory (which for the reasons stated above is
irrelevant to blood and is based on flawed and circular reasoning), the
authors describe a confusing ensemble of experimental measurements that do not
provide much valuable or relevant results and do not support their claims.

  * •

The measurements of the heat capacity ratio (HCR) of blood provided in Table 1
are completely unphysical, ranging from 5 to 118. No material even the most
exotic one exhibits such values and even for ideal gases, the equipartition
theorem states that the maximum value is $5/3\simeq 1.667$ for monoatomic gas.
Values shown in Table 1 are therefore wrong.

  * •

A confusing sentence in the methodology section says "We observed that around
60-85∘C, all the blood samples of human being boiloff in a non-linear
fashion". Blood, which is made mainly of water, does not boil at these
temperatures.

  * •

Similarly, in section II, they say "It is crystal clear from the case report
of Razavi et al. that the COVID-19 patient suffers from gas embolism because
the patient’s temperature exceeds 37.5 ∘ C". This entails two serious
problems, the first one being that the cited publication (reference 8 in their
publication list) does not make any mention of gas embolism in the reported
clinical case, it is therefore pure speculation from the authors. Then, they
seem to imply in this sentence (and at other instances in the paper) that
blood starts boiling or evaporating at temperatures above 37.5 ∘C, Which, had
it been true, would not have allowed the development of life on Earth as we
know it. Fortunately, the boiling point of blood (mainly composed of water) is
close to 100∘C at atmospheric pressure. Gas embolisms occur only in two
situations, either following accidental injection of gas bubbles in the
circulation (i.e. during surgery) or due to degassing of dissolved gases (N2,
CO2) following sudden and large decompression (i.e. underwater divers
ascending too quickly to the surface). The reported COVID-19 case is obviously
not concerned.

  * •

The _in vivo_ experiment of Fig. 12 is perplexing: they injected high-pressure
air into a rat’s artery to induce air embolism. It is hard to understand what
they were trying to prove: no quantitative pressure and velocity measurements
were made and the flow rate imposed by the syringe pump is not specified
either. Despite the authors’ claim, the numerical simulations of Fig. 10 and
11 are not correlated to the experiment at all since no parameter value
corresponds to experimental ones. This cannot constitute a proof of the claims
of the paper.

  * •

Figures 17, 18 and 19 only show compressed air being blown through a plastic
tube but do not demonstrate anything either qualitatively (no specific
observations are described) or quantitatively since no quantitative
measurements are made in these experiments. Besides, the connection between
these situations and blood flow in arteries is virtually inexistent. And there
is not even a demonstration of supersonic flow in these experiments.

  * •

More generally the focus on air embolisms in the second half of the paper
highlights a logical fallacy in the study: initially, the authors discuss
blood flow, attempt to measure its calorimetric properties in order to support
their theory, but in the end they claim to make experimental demonstrations
with air.

In conclusion, the experimental investigations performed by the authors and
comments based on other studies from the literature do not support the claims
derived from the model, which itself is based on wrong assumptions.

## References

  * Sanal Kumar  _et al._ (2022a) V. R. Sanal Kumar, B. R. Sundararam, P. K. Radhakrishnan, N. Chandrasekaran, S. K. Choudhary, V. Sankar, A. Sukumaran, V. Rajendran, S. A. R. M. Rafic, D. Panchal, Y. Raj, S. Shrivastava, C. Oommen, A. Jayaraman, D. Rajamanickam, and B. Srinivasan, “In vitro prediction of the lower/upper-critical biofluid flow choking index and in vivo demonstration of flow choking in the stenosis artery of the animal with air embolism,” Physics of Fluids 34, 101302 (2022a), https://doi.org/10.1063/5.0105407 .
  * Sanal Kumar  _et al._ (2022b) V. R. Sanal Kumar, V. Sankar, N. Chandrasekaran, V. Saravanan, A. Sukumaran, V. Rajendran, S. K. Choudhary, P. K. Radhakrishnan, C. Oommen, D. Panchal, S. A. R. M. Rafic, R. S. Bharath, Y. Raj, and S. Shrivastava, “Universal benchmark data of the three-dimensional boundary layer blockage and average friction coefficient for in silico code verification,” Physics of Fluids 34, 041301 (2022b), https://doi.org/10.1063/5.0086638 .
  * Kumar  _et al._ (2018) V. R. S. Kumar, V. Sankar, N. Chandrasekaran, V. Saravanan, V. Natarajan, S. Padmanabhan, A. Sukumaran, S. Mani, T. Rameshkumar, H. S. Nagaraju Doddi, K. Vysaprasad, S. Sharan, P. Murugesh, S. G. Shankar, M. N. Nejaamtheen, R. V. Baskaran, S. A. Rahman Mohamed Rafic, U. Harisrinivasan, and V. Srinivasan, “A closed-form analytical model for predicting 3d boundary layer displacement thickness for the validation of viscous flow solvers,” AIP Advances 8, 025315 (2018).
  * Fine and Millero (1973) R. A. Fine and F. J. Millero, “Compressibility of water as a function of temperature and pressure,” J. Chem. Phys. 59, 5529–5536 (1973), https://doi.org/10.1063/1.1679903 .
  * Robertson, Sequeira, and Kameneva (2008) A. M. Robertson, A. Sequeira, and M. V. Kameneva, “Hemorheology,” in  _Hemodynamical Flows: Modeling, Analysis and Simulation_ (Birkhäuser Basel, Basel, 2008) pp. 63–120.
  * Podgorski (2022) T. Podgorski, “Blood flow in large vessels,” (ISTE/WILEY, London/Hoboken, 2022) Chap. 1, pp. 1–32.
  * Duck (1990) F. A. Duck,  _Physical Properties of Tissues_ (Academic Press, London, 1990).
  * Lui  _et al._ (2020) M. Lui, S. Martino, M. Salerno, and M. Quadrio, “On the turbulence modeling of blood flow in a stenotic vessel,” J. Biomech. Eng.-Transactions Of The ASME 142 (2020), 10.1115/1.4044029.
  * Sanal Kumar  _et al._ (2021) V. R. Sanal Kumar, V. Saravanan, V. Srinivasan, S. Ganesh Shankar, S. Mani, V. Sankar, D. Krishnamoorthy, D. Natarajan, S. A. R. Mohamed Rafic, R. V. Baskaran, N. Chandrasekaran, V. Natarajan, A. Sukumaran, H. S. N. Doddi, S. Padmanabhan, G. Balasubramaniam, A. Saravanan, N. Santhosh Kumar, C. Lenin, A. Mohamed Imran Khan, M. Saravanan, R. Balakrishnan, R. Balachandru, M. Murugesan, D. V. Venkatesan, K. E. Shanjay, H. Sujith Kumar, N. A. Abhilash, D. Aswin Ram, V. Janardhanan, K. Krishnaraj, S. Gunasekaran, S. Karunakaran, S. Vineeshwar, A. Mariappan, R. Kiridharan, R. Sabarinathan, M. Ajith Kumar, R. Pradeep, U. K. Thianesh, M. Rajaram Perumal, K. S. Sumanth Eswar, and M. Rajesh, “The theoretical prediction of the boundary-layer-blockage and external flow choking at moving aircraft in ground effects,” Physics of Fluids 33, 036108 (2021), https://doi.org/10.1063/5.0040440 .
  * Ražnjević (1995) K. Ražnjević,  _Handbook of Thermodynamic Tables_ (Begell House, New York, 1995).
  * Mcintosh and Anderson (2010) R. L. Mcintosh and V. Anderson, “A comprehensive tissue properties database provided for the thermal assessment of a human at rest,” Biophysical Reviews and Letters 05, 129–151 (2010).