# Mixed Precision FGMRES-Based Iterative Refinement for Weighted Least Squares

Erin Carson111Department of Numerical Mathematics, Faculty of Mathematics and
Physics, Charles University, {carson<EMAIL_ADDRESS>Both authors
are supported by the Charles University GAUK project No. 202722, the Exascale
Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department
of Energy Office of Science and the National Nuclear Security Administration,
and by the European Union (ERC, inEXASCALE, 101075632). Views and opinions
expressed are those of the authors only and do not necessarily reflect those
of the European Union or the European Research Council. Neither the European
Union nor the granting authority can be held responsible for them. and Eda
Oktay11footnotemark: 1

#### Abstract

With the recent emergence of mixed precision hardware, there has been a
renewed interest in its use for solving numerical linear algebra problems fast
and accurately. The solution of least squares (LS) problems
$\min_{x}\|b-Ax\|_{2}$, where $A\in\mathbb{R}^{m\times n}$, arise in numerous
application areas. Overdetermined standard least squares problems can be
solved by using mixed precision within the iterative refinement method of
Björck, which transforms the least squares problem into an $(m+n)\times(m+n)$
”augmented” system. It has recently been shown that mixed precision GMRES-
based iterative refinement can also be used, in an approach termed GMRES-LSIR.
In practice, we often encounter types of least squares problems beyond
standard least squares, including weighted least squares (WLS),
$\min_{x}\|D^{1/2}(b-Ax)\|_{2}$, where $D^{1/2}$ is a diagonal matrix of
weights. In this paper, we discuss a mixed precision FGMRES-WLSIR algorithm
for solving WLS problems using two different preconditioners.

##  1 Introduction

Consider the weighted least squares problem

$\min_{x}\|D^{1/2}(Ax-b)\|_{2},$ (1)

where $A$ is an $m\times n$ matrix with $m\geq n$ and $D^{1/2}$ is a diagonal
matrix of weights. In applications where the weights vary significantly in
magnitude, transformation of the above into a standard least squares problem
is not numerically stable [2, Sec. 4.1.1]. Note that in some applications, the
weight matrix may be very ill-conditioned, for example, electrical networks,
certain finite element problems, and interior point methods [2, Remark 4.4.1].

Mathematically, the solution of (1) is given by the normal equations

$A^{T}DAx=A^{T}Db.$

Letting $y=D(b-Ax)$, this is also equivalent to the augmented system

$\begin{bmatrix}\alpha D^{-1}&A\\\
A^{T}&0\end{bmatrix}\begin{bmatrix}\alpha^{-1}y\\\
x\end{bmatrix}=\begin{bmatrix}b\\\ 0\end{bmatrix},$ (2)

where the scaling factor $\alpha$ has been introduced for stability; see [2,
Sec. 4.3.2 and Sec. 2.5.3].

One potential way to solve weighted problems is based on Householder QR
factorization, although one must use a Householder QR algorithm with row and
column pivoting in order to guarantee a good accuracy [2, Sec. 4.3.3]. Powell
and Reid investigated this in 1969 [8], and a reworked backward error analysis
in modern notation is given by Cox and Higham [5].

Another approach to solving such problems is using a mixed precision iterative
refinement (IR) algorithm. With the advent of modern GPUs which feature
multiple different hardware precisions, this approach has seen renewed
interest. See Table 1 for various IEEE precisions and their units roundoff.
The general scheme of IR introduced in [13] is given in Algorithm 1. Depending
on the precision $u_{f}$ chosen for factorization $u_{r}$ for the residua
computation, $u_{s}$ for the correction solve, and the working precision $u$
for storing data and the solution, one gets a variant of the mixed precision
IR algorithm. To solve LS problems, the author in [1] used a 2-precision IR
approach to solve the augmented system; we call this LS-IR. To reduce the
computation and memory cost, LS-IR never forms the augmented system
explicitly, using only the QR factors of $A$. However, the QR factorization
can be very expensive if $A$ is very large. Thus, the authors in [4] use three
precisions in LS-IR to further reduce the computation cost. Using a lower
precision $u_{f}$ in the QR factorization can provide a cheaper algorithm
while maintaining accuracy. Based on the analysis in [4], we can say as long
as $\kappa_{\infty}(\widetilde{A})\lesssim u_{f}^{-1}$, the backward error of
three precision LS-IR is $\mathcal{O}(u)$ and the forward error is

$\dfrac{\|\tilde{x}-\hat{\tilde{x}}\|_{\infty}}{\|\tilde{x}\|_{\infty}}\approx
u_{r}cond(\tilde{A},\tilde{x})+u.$ (3)

Above, (3) shows that if $u_{f}$ is chosen to be fp16, LS-IR converges only
for $\kappa_{\infty}(\tilde{A})\lesssim 4.88\cdot 10^{4}$, which restricts the
set of problems for which the algorithm is guaranteed to converge.

Algorithm 1 IR [13]

1:$A\in\mathbb{R}^{m\times n}$, $(m\geq n)$, $b\in\mathbb{R}^{m}$

2:$Ax_{0}=b$ in precision $u_{f}$; store in precision $u$

3:for $i=0$: $i_{max}-1$ do

4: $r_{i}=b-Ax_{i}$ in precision $u_{r}$; store in precision $u$

5: $Ad_{i+1}=r_{i}$ in precision $u_{s}$; store in precision $u$

6: $x_{i+1}=x_{i}+d_{i+1}$ in precision $u$

7: if converged then return $x_{i+1}$

Table 1: Various IEEE precisions and their units roundoff.  Precision | Unit Roundoff
---|---
fp16 (half) | $4.88\cdot 10^{-4}$
fp32 (single) | $5.96\cdot 10^{-8}$
fp64 (double) | $1.11\cdot 10^{-16}$
fp128 (quad) | $9.63\cdot 10^{-35}$

##  2 GMRES-based approach

To allow the use of low precision factorization for more ill-conditioned
systems, we can adapt the GMRES-based approach of [4] for non-weighted least
squares problems, called GMRES-LSIR. The key difference with LS-IR is that
this approach uses preconditioned GMRES to solve the linear system in line 5
of Algorithm 1.

For the non-weighted case, i.e., $D=I$, the authors in [4] propose the left
preconditioner

$M=\begin{bmatrix}\alpha I&Q\widehat{R}\\\
\widehat{R}^{T}Q^{T}&0\end{bmatrix},$ (4)

composed of the QR factors of $A$ denoted as $Q$ and $\widehat{R}$.

According to [4], as long as $\kappa_{\infty}(A)\leq u^{-1/2}u_{f}^{-1}$,
GMRES-LSIR has $\mathcal{O}(u)$ backward and forward error. In the simplified
case where $u=u_{r}$, the forward error is (3) if $\kappa_{\infty}(A)\leq
u^{-1/3}u_{f}^{-2/3}$. Using the preconditioner $M$, the condition number of
the preconditioned augmented system $M^{-1}\tilde{A}$ can be bounded by

$\kappa_{\infty}(M^{-1}\tilde{A})\lesssim(1+2m\sqrt{n}\tilde{\gamma}_{mn}^{f}\kappa_{\infty}(A))^{2},\quad\text{with
}\quad\tilde{\gamma}_{mn}^{f}=\dfrac{cmn}{1-mnu_{f}}$ (5)

for a small constant $c$. This bound shows that even if $\kappa_{\infty}(A)\gg
u_{f}^{-1}$, $M$ computed in precision $u_{f}$ will reduce
$\kappa_{\infty}(M^{-1}\tilde{A})$, which is one of the most crucial benefits
of GMRES-LSIR, allowing more ill-conditioned problems to be solved.

To solve WLSP using the iterative refinement-based approach, we can generally
follow the same approach as described above. As in [4, Sec. 3.1], we want to
develop a preconditioner $M$ for our approach. Our aim is thus to find an
effective and inexpensive preconditioner $M$ and prove that the resulting
preconditioned coefficient matrix is sufficiently well-conditioned to
guarantee the convergence of iterative refinement. In this study, we focus on
two different preconditioners: a left preconditioner, which is a direct
extension of the above, and a split preconditioner.

##  3 Left QR Preconditioner

Consider the preconditioner

$M_{l}=\begin{bmatrix}\alpha D^{-1}&Q\widehat{R}\\\
\widehat{R}^{T}Q^{T}&0\end{bmatrix},$ (6)

which is the direct extension of the $M$ in (4), except with a $D^{-1}$
instead of $I$ in the $(1,1)$ block. We can explicitly write down the inverse
of this matrix as

$M_{l}^{-1}=\begin{bmatrix}\frac{1}{\alpha}\left(D-DQ(Q^{T}DQ)^{-1}Q^{T}D\right)&DQ(Q^{T}DQ)^{-1}\widehat{R}^{-T}\\\
\widehat{R}^{-1}(Q^{T}DQ)^{-1}Q^{T}D&-\alpha\widehat{R}^{-1}(Q^{T}DQ)^{-1}\widehat{R}^{-T}\end{bmatrix}.$

Let

$\widetilde{E}=\widetilde{A}-M_{l}=\begin{bmatrix}0&-E\\\
-E^{T}&0\end{bmatrix},$

where the error term $E$ is defined as

$A+E=Q\widehat{R},$ (7)

due to finite precision computation of the QR factorization in some precision
$u_{f}$ [5].

Now, note that we can write

$\displaystyle M_{l}^{-1}\widetilde{A}$
$\displaystyle=M_{l}^{-1}(M_{l}+\widetilde{E})=I+M_{l}^{-1}\widetilde{E}$
$\displaystyle\widetilde{A}^{-1}M_{l}$
$\displaystyle=(M_{l}+\widetilde{E})^{-1}M_{l}\approx
I-M_{l}^{-1}\widetilde{E}.$

Thus

$\kappa_{\infty}(M_{l}^{-1}\widetilde{A})=\|M_{l}^{-1}\widetilde{A}\|_{\infty}\|\widetilde{A}^{-1}M_{l}\|_{\infty}\lesssim(1+\|M_{l}^{-1}\widetilde{E}\|_{\infty})^{2},$
(8)

so if we bound the quantity $\|M_{l}^{-1}\widetilde{E}\|_{\infty}$, we obtain
a bound on the condition number of the preconditioned system. We can
explicitly write

$M_{l}^{-1}\widetilde{E}=\begin{bmatrix}-DQ(Q^{T}DQ)^{-1}\widehat{R}^{-T}E^{T}&-\frac{1}{\alpha}\left(D-DQ(Q^{T}DQ)^{-1}Q^{T}D\right)E\\\
\alpha\widehat{R}^{-1}(Q^{T}DQ)^{-1}\widehat{R}^{-T}E^{T}&-\widehat{R}^{-1}(Q^{T}DQ)^{-1}Q^{T}DE\end{bmatrix}.$

Note that if $D=I$ this reduces to the case in [4] and gives (4).

The quantity $(QDQ)^{-1}Q^{T}D$ is what is called the “scaled” or “weighted”
pseudoinverse [11]; see also [9], [12]. In essence, the results of Stewart
[11] show that this quantity is bounded and is independent of $D$. We can,
however, still give a concrete bound on this quantity, which may in practice
be a significant overestimate. Let $Q_{D}=D^{1/2}Q$. Then

$\displaystyle\|(Q^{T}DQ)^{-1}Q^{T}D\|_{2}$
$\displaystyle=\|(Q_{D}^{T}Q_{D})^{-1}Q_{D}^{T}D^{1/2}\|_{2}$ (9)
$\displaystyle\leq\|Q_{D}^{\dagger}\|_{2}\|D^{1/2}\|_{2}$
$\displaystyle\leq\|Q^{\dagger}\|_{2}\|D^{-1/2}\|_{2}\|D^{1/2}\|_{2}$
$\displaystyle=\kappa_{2}^{1/2}(D).$

As an alternative, we can let

$\|(Q^{T}DQ)^{-1}Q^{T}D\|_{2}\equiv\rho,$

and then argue that this $\rho$ is hopefully of reasonable size since it is
independent of $D$, using the results of Stewart (we will do this for now).

Now, we want to derive a bound for $\|M_{l}^{-1}\widetilde{E}\|_{\infty}$.
Following [4], we have

$\displaystyle\|M_{l}^{-1}\widetilde{E}\|_{\infty}$
$\displaystyle\leq\max(\|(M_{l}^{-1}\widetilde{E})(1,1)\|_{\infty}+\|(M_{l}^{-1}\widetilde{E})(1,2)\|_{\infty},$
$\displaystyle\phantom{\leq}\qquad\qquad\qquad,\|(M_{l}^{-1}\widetilde{E})(2,1)\|_{\infty}+\|(M_{l}^{-1}\widetilde{E})(2,2)\|_{\infty})$
$\displaystyle\leq\sqrt{m}\max\left(\|(M_{l}^{-1}\widetilde{E})(1,1)\|_{F},\|(M_{l}^{-1}\widetilde{E})(2,1)\|_{F}\right)$
$\displaystyle\phantom{\leq}+\sqrt{n}\max\left(\|(M_{l}^{-1}\widetilde{E})(1,2)\|_{F},\|(M_{l}^{-1}\widetilde{E})(2,2)\|_{F}\right),$
(10)

where $(M_{l}^{-1}\widetilde{E})(i,j)$ denotes the $(i,j)$-block of
$M_{l}^{-1}\widetilde{E}$.

Using $|\alpha|\approx\|A^{\dagger}\|_{2}^{-1}$ and ignoring terms of order
$u_{f}^{2}$,

$\displaystyle\|(M_{l}^{-1}\widetilde{E})(1,1)\|_{F}$
$\displaystyle=\|DQ(Q^{T}DQ)^{-1}\widehat{R}^{-T}E^{T}\|_{F}\leq\rho\|A^{\dagger}\|_{2}\|E^{T}\|_{F}$
$\displaystyle\|(M_{l}^{-1}\widetilde{E})(1,2)\|_{F}$
$\displaystyle=\left\|\frac{1}{\alpha}D\left(I-Q(Q^{T}DQ)^{-1}Q^{T}D\right)E\right\|_{F}\leq\rho\|D\|_{2}\|A^{\dagger}\|_{2}\|E\|_{F}$
$\displaystyle\|(M_{l}^{-1}\widetilde{E})(2,2)\|_{F}$
$\displaystyle=\|\widehat{R}^{-1}(Q^{T}DQ)^{-1}Q^{T}DE\|_{F}\leq\rho\|A^{\dagger}\|_{2}\|E\|_{F}.$

For the final term, we have

$\displaystyle\|(M_{l}^{-1}\widetilde{E})(2,1)\|_{F}$
$\displaystyle=\|\alpha\widehat{R}^{-1}(Q^{T}DQ)^{-1}\widehat{R}^{-T}E^{T}\|_{F}$
$\displaystyle=\|\alpha\widehat{R}^{-1}(Q^{T}DQ)^{-1}Q^{T}D(Q^{T}D)^{\dagger}\widehat{R}^{-T}E^{T}\|_{F}$
$\displaystyle\leq\rho\|D^{-1}\|_{2}\|A^{\dagger}\|_{2}\|E^{T}\|_{F}.$

Then plugging into (10), and letting
$\theta\equiv\max(1,\|D\|_{2},\|D^{-1}\|_{2})$, we have

$\displaystyle\|M_{l}^{-1}\widetilde{E}\|_{\infty}$
$\displaystyle\leq(\sqrt{m}+\sqrt{n})\rho\theta\|A^{\dagger}\|_{2}\|E\|_{F}.$

Assuming that we have a standard QR factorization,
$\|E\|_{F}\leq\tilde{\gamma}_{mn}^{f}\|A\|_{F}$, and thus

$\displaystyle\|M_{l}^{-1}\widetilde{E}\|_{\infty}$
$\displaystyle\leq(\sqrt{m}+\sqrt{n})\rho\theta\|A^{\dagger}\|_{2}\|E\|_{F}$
(11) $\displaystyle\leq
2m\sqrt{n}\rho\theta\tilde{\gamma}_{mn}^{f}\kappa_{\infty}(A).$ (12)

Note that in the case that we use a standard QR factorization and $D=I$ (and
thus $\rho=1$), we recover the same bound as in [4, Section 3.1]. Finally,
plugging into (8), we obtain

$\kappa_{\infty}(M_{l}^{-1}\widetilde{A})\lesssim\left(1+2m\sqrt{n}\rho\theta\tilde{\gamma}_{mn}^{f}\kappa_{\infty}(A)\right)^{2}.$
(13)

##  4 Block-diagonal Split Preconditioner

Note that the dependence of the bound in (11) on $D$ is not ideal since $D$
can be very ill-conditioned. For the behavior of the condition number of
$M_{l}^{-1}\widetilde{A}$, see Section 5.1. To construct a preconditioner
which results in a preconditioned matrix whose conditioning does not depend on
$D$, we thus consider using the block diagonal preconditioner

$M_{b}=\begin{bmatrix}\alpha D^{-1}&0\\\ 0&\hat{C}\end{bmatrix},$ (14)

where $\hat{C}$ is a symmetric positive definite approximation to the Schur
complement, $\alpha^{-1}A^{T}DA$ [10].

Although $M_{b}$ can be used as a left/right preconditioner, to obtain a
symmetric preconditioned system, we consider using it as a split
preconditioner and study the system

$\displaystyle M_{b}^{-1/2}\tilde{A}M_{b}^{-1/2}$
$\displaystyle=\begin{bmatrix}(\alpha D^{-1})^{-1/2}&0\\\
0&\hat{C}^{-1/2}\end{bmatrix}\begin{bmatrix}\alpha D^{-1}&A\\\
A^{T}&0\end{bmatrix}\begin{bmatrix}(\alpha D^{-1})^{-1/2}&0\\\
0&\hat{C}^{-1/2}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}(\alpha
D^{-1})^{-1/2}(\alpha D^{-1})(\alpha D^{-1})^{-1/2}&(\alpha
D^{-1})^{-1/2}A^{T}\hat{C}^{-1/2}\\\ \hat{C}^{-1/2}A^{T}(\alpha
D^{-1})^{-1/2}&0\end{bmatrix}$ $\displaystyle=\begin{bmatrix}I&\hat{A}\\\
\hat{A}^{T}&0\end{bmatrix}.$

No applicable analysis of the forward and backward errors in split
preconditioned GMRES exists, but there is an existing analysis of these
quantites for a general split preconditioned flexible GMRES method (FGMRES).
Thus, to be able to use split preconditioners and discuss their effects on the
stability of this general approach, we use FGMRES instead of GMRES in GMRES-
LSIR. Our new variant to solve WLSP is therefore called FGMRES-WLSIR.

We use a symmetric positive definite approximation $\hat{C}$ such that it is
spectrally equivalent to the matrix $\alpha^{-1}A^{T}DA$, i.e., there exist
positive constants $0<\hat{\gamma}\leq\hat{\delta}$ such that

$\hat{\gamma}(\hat{C}y,y)\leq(B^{T}(\alpha
D^{-1})^{-1}Ay,y)\leq\hat{\delta}(\hat{C}y,y)$

for all vectors $y\in\mathbb{R}^{n}$. Using [10, Proposition 3.4], we can find
the spectrum of the preconditioned matrix:

$sp(M_{b}^{-1/2}\tilde{A}M_{b}^{-1/2})=\frac{1}{2}\left(1\pm\sqrt{1+4\sigma^{2}_{k}(\hat{A})}\right)\quad
for\quad k=1,\ldots,n-r,$

where rank($\hat{A}$) = $n-r$ and
$0<\sigma_{n-r}(\hat{A})\leq\cdots\leq\sigma_{1}(\hat{A})$ are the singular
values of $\hat{A}$.

We assume that the Schur complement is computed exactly and we consider the
left-preconditioned matrix

$M_{b}^{-1}\tilde{A}=\begin{bmatrix}I&\alpha^{-1}A^{T}D\\\
\alpha^{-1}A(A^{T}DA)^{-1}&0\end{bmatrix},$

which is a nonsymmetric diagonalizable matrix with three distinct eigenvalues
$\left\\{1,\frac{1}{2}(1\pm\sqrt{5})\right\\}$. This makes the block diagonal
preconditioner a suitable choice for FGMRES since the method can converge in a
small number of iterations. However, our experiments show that
$\kappa_{\infty}(M_{b}^{-1}\tilde{A})$ can be very large, often larger than
$\kappa_{\infty}(\tilde{A})$, when $A$ is ill-conditioned. An ill-conditioned
preconditioned matrix makes the preconditioner unsuitable for proving the
backward stability of FGMRES although we show in Section 5 that in practice,
the split preconditioner improves the condition number of the preconditioned
matrix.

Using the analysis in [10], we can obtain the bound

$\kappa_{\infty}(M_{b}^{-1/2}\tilde{A}M_{b}^{-1/2})\leq\dfrac{\hskip
1.0pt|1+\sqrt{1+4\sigma_{1}^{2}(\hat{A})}|\hskip 1.0pt}{\hskip
1.0pt|1-\sqrt{1+4\sigma_{n-r}^{2}(\hat{A})}|\hskip 1.0pt}(n+m),$ (15)

where $0<\sigma_{n-r}(\hat{A})\leq\ldots\leq\sigma_{1}(\hat{A}),$

$\hat{A}=\dfrac{1}{\sqrt{\alpha}}D^{1/2}A^{T}\hat{C}^{-1/2}\text{ with
}\hat{C}\approx\alpha^{-1}A^{T}DA.$

Although there is still a dependence on $D$ in (15), we can numerically
eliminate its impact by computing the Schur complement in a specific manner.
For $\hat{C}$, we compute $D^{1/2}A$ in high precision and use the R-factor
(from lower precision QR factorization) of it in $M_{b}^{1/2}$. The numerical
experiments in Section 5 show that this trick reduced the effect of the
conditioning of $D$.

On the other hand, using a split preconditioner in (F)GMRES has several
disadvantages regarding error analysis. Using the analysis in [3], we see that
for the block diagonal split preconditioner (the left and right
preconditioners are both $M_{b}^{1/2}$), the forward error of FGMRES is
bounded by

$\dfrac{\|x-\bar{x}_{k}\|}{\|x\|}\lesssim\kappa_{\infty}(M_{b}^{-1/2}\tilde{A}M_{b}^{-1/2})\kappa_{\infty}(M_{b}^{1/2})\mathcal{O}(u),$
(16)

whereas for any left preconditioner, the bound is given as

$\dfrac{\|x-\bar{x}_{k}\|}{\|x\|}\lesssim\kappa_{\infty}(M_{l}^{-1}\tilde{A})\mathcal{O}(u).$

Even though $M_{b}^{-1/2}\tilde{A}M_{b}^{-1/2}$ is well-conditioned,
$M_{b}^{1/2}$ can still be ill-conditioned, which can cause the forward error
to be more than $1$ and thus cannot guarantee FGMRES-WLSIR convergence. This
result shows the need for a bound on $\kappa_{\infty}(M_{b}^{1/2})$.

##  5 Numerical Experiments

To illustrate the impact of different preconditioners on FGMRES-WLSIR
convergence, we perform several numerical experiments in MATLAB. The
experiments are performed on a computer with AMD Ryzen 5 4500U having 6 CPUs
and 8 GB RAM with OS system Ubuntu 22.04 LTS. In our numerical experiments, we
used built-in MATLAB datatypes for double and single precisions. To simulate
half-precision floating point arithmetic, we use the chop library and
associated functions from [7] available at `https://github.com/higham/chop`
and `https://github.com/SrikaraPranesh/`
`LowPrecision_Simulation`. Code for our FGMRES-WLSIR and associated functions
can be found in the repository `https://github.com/edoktay/fgmreswlsir`.

For the examples in this section, we set $A^{\prime}=$ gallery(‘randsvd’,
[100,10], 1e2, 3), which creates a matrix with 2-norm condition number
$10^{2}$ of size $100\times 10$ with geometrically distributed singular
values. We prescale the rows of $A^{\prime}$ so that they have drastically
different sizes. We use 9 different scalings $A=S^{-1}A^{\prime}$ where $S$ is
a diagonal matrix created by diag(logspace(1, j, 100) where
$j\in[1,2,4,6,8,10,12,14,16]$. We set the weight matrix $D$ such that each row
$i$ of $A$ has $\max|A(i,j)|=1$. The scaling factor
$\alpha=2^{-1/2}\sigma_{n}$, where $\sigma_{n}$ is the smallest singular value
of $A$. We compute QR factorizations in half, single, and double precision and
construct the corresponding preconditioners $M$. We then measure the infinity-
norm condition number of the preconditioned system.

In each figure, the condition number of the preconditioned systems is
represented as colored lines. Each color shows half (red), single (green), and
double (blue) precisions used for computing QR factorizations for the
preconditioners. The dashed black line gives the condition number of the
unpreconditioned augmented system and the dotted black line gives the inverse
of the unit roundoff for the FGMRES-WLSIR working precision $u$. The figure
shows that the convergence of FGMRES-WLSIR is guaranteed only when the colored
solid lines remain below the dashed line. Numerically, this shows the cases
when $\kappa_{\infty}(M^{-1}\widetilde{A})\leq u^{-1}$. Only in this case can
we guarantee that the forward error of FGMRES is less than $1$ and thus
FGMRES-WLSIR converges.

###  5.1 Dependence of $\|M_{l}^{-1}\widetilde{E}\|_{\infty}$ on $D$

We perform a numerical experiment to demonstrate how
$\kappa_{\infty}(M_{l}^{-1}\widetilde{A})$ changes with the conditioning of
$D$. For this example, we assume that FGMRES-WLSIR is performed using single
precision as the working precision $u$. The results are shown in Figure 1.

Figure 1: Measured condition number of the preconditioned systems and
estimates of the bound on the condition number (13) where the preconditioners
(6) are constructing using QR factorizations computed in various precisions,
versus the condition number of the weight matrix $D$. The working precision
for FGMRES-WLSIR is assumed to be single precision.

The figure shows that using $M_{l}$ in double precision preserves stability
even if $D$ is very badly-conditioned. When we switch to single precision, we
see $M_{l}$ is useful (convergence of FGMRES is guaranteed) when
$\kappa_{\infty}(D)<10^{10}$. On the other hand, using half-precision limits
the usability of the preconditioner since the condition number of $D$ should
be at most $10^{6}$, which may not be the case in real-world problems. We thus
conclude that using lower precision can be tricky for $M_{l}$ due to the fact
that its error bound depends on the conditioning of $D$. We also observe from
the black dashed line that half-precision $M_{l}$ does not change the
conditioning of the preconditioned system significantly.

###  5.2 $M_{l}$ versus $M_{b}$

To examine the effect of preconditioners on the conditioning of the augmented
system, we construct Table 2 using random dense matrices with a randomly
generated solution vector $b$. We construct matrix $A$ in a similar manner as
in Section 5.1. The table shows that the preconditioner $M_{l}$ does not
decrease the condition number sufficiently, whereas $M_{b}$ works well. We
finally observe from the last column that even though $M_{b}^{1/2}$ is ill-
conditioned since the preconditioned system is very well conditioned, FGMRES-
WLSIR still converges due to the forward error constraint in (16).

Table 2: Condition numbers of $\widetilde{A}$, right preconditioner, and preconditioned augmented matrices. $\kappa_{2}(A)$ | $\kappa_{\infty}(\widetilde{A})$ | $\kappa_{\infty}(M_{l}^{-1}\widetilde{A})$ | $\kappa_{\infty}(M_{b}^{-1/2}\widetilde{A}M_{b}^{-1/2})$ | $\kappa_{\infty}(M_{b}^{1/2})$
---|---|---|---|---
1.00e+02 | 1.12e+04 | 3.73e+00 | 8.56e+01 | 1.05e+03
1.12e+02 | 4.92e+03 | 2.55e+00 | 7.80e+01 | 3.96e+02
1.47e+02 | 1.73e+05 | 1.16e+02 | 5.79e+01 | 3.60e+02
1.91e+02 | 8.04e+08 | 1.14e+05 | 4.03e+01 | 3.08e+04
2.47e+02 | 5.06e+12 | 1.50e+08 | 3.02e+01 | 2.48e+06
3.21e+02 | 3.81e+16 | 4.36e+11 | 2.62e+01 | 2.16e+08
4.22e+02 | 2.97e+20 | 5.24e+14 | 2.30e+01 | 1.87e+10
5.61e+02 | 2.26e+24 | 1.12e+18 | 2.05e+01 | 1.63e+12
7.59e+02 | 1.67e+28 | 1.16e+22 | 1.86e+01 | 1.36e+14

To demonstrate the difference in the numerical behavior of the preconditioned
system with $M_{l}$ and $M_{b}$ with the conditioning of $D$, we also perform
several numerical experiments using `ash958` and `robot24c1_mat5` matrices
from the SuiteSparse collection [6]. We prescale the rows of $A^{\prime}$ so
that they have drastically different sizes. We use 9 different scalings
$A=SA^{\prime}$ where $S$ is the scaling matrix used in Section 5.1. The
properties of both matrices are given in Table 3. For the experiments in this
section, we assume that we use double precision as the working precision $u$
in FGMRES-WLSIR.

Table 3: Properties of matrices from the SuiteSparse collection.  Name | $m$ | $n$ | $\kappa_{2}(A)$ |  $\\#$nnz
---|---|---|---|---
`ash958` | 958 | 292 | 3.2014 | 1916
`robot24c1_mat5` | 404 | 302 | $3.33\times 10^{11}$ | 15118

Figure 2: Measured condition number of the preconditioners (left) and
preconditioned systems (right) using ash958 matrix as $A$ where $M_{l}$ and
$M_{b}$ are constructed using QR factorizations in various precisions, versus
the condition number of the weight matrix $D$. The working precision for
FGMRES-WLSIR is assumed to be double precision.

We observe from the right plot in Figure 2 that single precision $M_{l}$ can
handle up to $\kappa_{\infty}(D)<10^{13}$ however, we note that the algorithm
can still be expensive due to QR factorization. On the other hand, we see from
the same plot that using $M_{b}$, even with a very ill-conditioned $D$ will
provide a well-conditioned preconditioned coefficient matrix. Numerically, we
don’t see the effect of $D$ because of our way of computing the Schur
complement. However, because of the dependence of the forward error of FGMRES
on $\kappa_{\infty}(M_{b}^{1/2})$, we need to look at the left plot as well.
From the left plot we see that even if numerically, the preconditioned system
is very well-conditioned, because of the conditioning of the right split
preconditioner (in our case, it is equivalent to the left split
preconditioner), $\kappa_{\infty}(M_{b}^{1/2})$ is sufficiently small only
when $\kappa_{\infty}(D)<10^{9}$. For `ash958`, we can say that using the
split block diagonal preconditioner $M_{b}$ does not give any significant
advantage over $M_{l}$ in fp16.

Figure 3: Measured condition number of the preconditioners (left) and
preconditioned systems (right) using robot24c1_mat5 matrix as $A$ where
$M_{l}$ and $M_{b}$ are constructed using QR factorizations in various
precisions, versus the condition number of the weight matrix $D$. The working
precision for FGMRES-WLSIR is assumed to be double precision.

From the right plot in Figure 3, we see that using half-precision in computing
$M_{l}$ is not useful even if $D$ is very well-conditioned. Moreover, the left
preconditioner works only when $\kappa_{\infty}(D)<10^{7}$. For the split
preconditioner, we again observe the same behavior as the previous example.
From the left preconditioner, we again expect a limitation coming from the
conditioning of $M_{b}^{1/2}$. We observe that in any precision,
$\kappa_{\infty}(M_{b}^{1/2})\leq u^{-1}$ when $\kappa_{\infty}(D)<10^{7}$.
However, even though $\kappa_{\infty}(M_{b}^{1/2})\leq u^{-1}$, since the
forward error bound is obtained via the multiplication of both condition
numbers, it will be close to 1 for $M_{b}$ in fp16 and fp32 when $D$ is well-
conditioned. Therefore, using fp16 in computing both preconditioners is not
applicable to this matrix. Furthermore, we see that for `robot24c1_mat5`,
$M_{l}$ is more useful than $M_{b}$ in terms of fp32 and fp64 applicability.

##  6 Conclusion

In various areas, problems may need a very badly-conditioned weight matrix to
be able to be modeled as least squares problems. Most of the available fast
least squares solvers are not directly usable for weighted cases. One thus
needs to make some changes in algorithms such as changing the preconditioning.
GMRES-LSIR is one of the potential iterative solvers in the literature for
solving least squares problems fast and accurately using lower precision. With
this motivation, our goal was to extend GMRES-LSIR to solve weighted least
squares problems. For this extension, we examined the use of two different
preconditioners, namely left and block split preconditioners, in the
algorithm. To construct an error bound for the split preconditioner, we use
FGMRES instead of the GMRES algorithm and introduce our approach, FGMRES-
WLSIR.

Using the analyses in the literature, we introduce error bounds for both
preconditioners under assumptions on the conditioning of $A$. From our
analysis, we observe that the dependence of the weight matrix in the error
bound of the left preconditioner limits the use of low precisions in our
approach. To reduce this dependence numerically, we study the block split
preconditioner. By computing the Schur complement in a specific way, we
numerically reduce the dependence on the conditioning of $D$. However, from
the preconditioned FGMRES analysis, we observe that the forward error of
FGMRES also depends on the conditioning of the right preconditioner in the
split preconditioner case. We therefore examine the conditioning of both
preconditioners and their ranges of applicability numerically. We conclude
that since the conditioning of the weight matrix is highly problem-dependent,
we cannot generalize which preconditioner is more useful for FGMRES-WLSIR.
Although our numerical experiments show that using the block split
preconditioner in half-precision may be used in more ill-conditioned systems,
in real applications, $D$ might be worse-conditioned.

Because of the dependence of both preconditioners on $D$ and the dependence of
split preconditioned FGMRES on the right preconditioner, further study is
warranted. Future studies can focus either on the choice of a preconditioner
or another iterative approach other than FGMRES. In any case, the optimal
approach will ideally have error bounds independent of the weight matrix.

## References

  * [1] Å. Björck, Iterative refinement of linear least squares solutions i, BIT Numerical Mathematics, 7 (1967), pp. 257–278.
  * [2] Å. Björck, Numerical methods for least squares problems, SIAM, 1996.
  * [3] E. Carson and I. Daužickaitė, The stability of split-preconditioned fgmres in four precisions, arXiv preprint arXiv:2303.11901, (2023).
  * [4] E. Carson, N. J. Higham, and S. Pranesh, Three-precision GMRES-based iterative refinement for least squares problems, SIAM Journal on Scientific Computing, 42 (2020), pp. A4063–A4083.
  * [5] A. J. Cox and N. J. Higham, Stability of Householder QR factorization for weighted least squares problems, in Numerical Analysis 1997, Proceedings of the 17th Dundee Biennial Conference, D. F. Griffiths, D. J. Higham, and G. A. Watson, eds., vol. 380 of Pitman Research Notes in Mathematics, Addison Wesley Longman, Harlow, Essex, UK, 1998, pp. 57–73.
  * [6] T. A. Davis and Y. Hu, The university of florida sparse matrix collection, ACM Trans. Math. Softw., 38 (2011), https://doi.org/10.1145/2049662.2049663, https://doi.org/10.1145/2049662.2049663.
  * [7] N. J. Higham and S. Pranesh, Simulating low precision floating-point arithmetic, SIAM Journal on Scientific Computing, 41 (2019), pp. C585–C602, https://doi.org/10.1137/19M1251308, https://doi.org/10.1137/19M1251308, https://arxiv.org/abs/https://doi.org/10.1137/19M1251308.
  * [8] A. Morell, ed., On applying Householder’s method to linear least squares problems, North-Holland, Amsterdam, 1969.  pp. 122-126.
  * [9] D. P. O’Leary, On bounds for scaled projections and pseudoinverses, Linear Algebra and its Applications, 132 (1990), pp. 115–117.
  * [10] M. Rozložník, Saddle-point problems and their iterative solution, Springer, 2018.
  * [11] G. W. Stewart, On scaled projections and pseudoinverses, Linear Algebra and its Applications, 112 (1989), pp. 189–193.
  * [12] M. Wei, Upper bound and stability of scaled pseudoinverses, Numerische Mathematik, 72 (1995), pp. 285–293.
  * [13] J. H. Wilkinson, Rounding errors in algebraic processes, Prentice-Hall, 1963.


11institutetext: Arizona State University, Tempe, AZ, USA 22institutetext:
University of Illinois Chicago, Chicago, Illinois, USA
22email<EMAIL_ADDRESS>
# Bridging the Gap: Commonality and Differences between Online and Offline
COVID-19 Data

Nayoung Kim 11 Ahmadreza Mosallanezhad 11 Lu Cheng 1 and 21 and 2 Baoxin Li 11
Huan Liu 11

###### Abstract

With the onset of the COVID-19 pandemic, news outlets and social media have
become central tools for disseminating and consuming information. Because of
their ease of access, users seek COVID-19-related information from online
social media (i.e., online news) and news outlets (i.e., offline news). Online
and offline news are often connected, sharing common topics while each has
unique, different topics. A gap between these two news sources can lead to
misinformation propagation. For instance, according to the Guardian, most
COVID-19 misinformation comes from users on social media. Without fact-
checking social media news, misinformation can lead to health threats. In this
paper, we focus on the novel problem of bridging the gap between online and
offline data by monitoring their common and distinct topics generated over
time. We employ Twitter (online) and local news (offline) data for a time span
of two years. Using online matrix factorization, we analyze and study online
and offline COVID-19-related data differences and commonalities. We design
experiments to show how online and offline data are linked together and what
trends they follow.

###### Keywords:

Matrix Factorization COVID-19 Social Media Analysis

##  1 Introduction

The COVID-19 pandemic has been accompanied by a massive infodemic involving a
variety of topics, such as vaccination, inflation, and so on. Communication of
important information during emergency situations is critical [17, 2] for
taking actions to contain disease. With its convenience, easy access, and
large volume, news media becomes an expected means to a pandemic response.
People seek and receive information from numerous sources such as newspapers,
television, and increasingly, social media. For example, during the early
stages of the COVID-19 pandemic, Twitter and news outlets such as Washington
post reported about US states’ social distancing mandates to help prevent the
growth of the virus.

Due to the different nature of online social media data and offline news
articles, information gaps may exist among online social media users and
offline news media readers. These users and readers consist of scientists,
medical and public health professionals, and the general public. A large
information gap could contribute to the spread of misinformation in news
media. Therefore, to better understand the relations between online social
media data and offline news articles, we aim to monitor their common and
distinct topics generated over time. In particular, connecting offline with
online data can be beneficial in two main ways: (1) Researchers can measure
public interest over time and identify new themes of COVID-19-related
discussion emerging quickly through data streams. (2) By providing users with
useful information about COVID-19, it reduces the chance of exposing them to
misinformation, i.e., a gap between offline and online topics can alert users
about possible misinformation propagation.

Figure 1 shows the problem setting. Our goal is to discover and link topics
from the traditional and social media news sources, to understand their
association. We study this problem on a COVID-19-related dataset [7] collected
through Twitter and Google News. This problem presents the following
challenges: (1) The collected dataset contains a significant amount of noisy
and irrelevant data. In the online data, not all the tweets follow a specific
topic or provide useful information about an actual event. Consequently, the
integration approach should discard noise and irrelevant data. (2) Identifying
topics across two different data sources in a separate and joint setting. The
model should be able to find the association between them by leveraging the
dataset properties. Typically, Twitter data are short texts, including
hashtags, user mentions, and covering any topics. In contrast, news articles
published on media websites are long and structured.

Figure 1: Illustration of the common and distinct topics identified within two
types of datasets: offline and online.

To address the challenges above, we use a variant of the online Non-negative
matrix factorization (NMF) method, joint ONMF. NMF methods are capable of
dimensionality reduction and data clustering simultaneously. Because of our
data’s nature and its streaming property, we cannot use vanilla NMF. Further,
the online NMF method cannot model discovered topics over time. As a result,
we adopt joint ONMF to find common and distinct topics over time and link
online and offline data. Specifically, we overcome the challenges by making
the following contributions: (1) We study the novel problem of linking data
related to COVID-19 from online social media and offline news sources, (2) We
use joint ONMF to identify common and distinct topics across time, and (3) We
conduct extensive experiments on real-world datasets to draw conclusions.

##  2 Related Work

There have been a growing number of research on social media discourse
associated with disasters or crises. For example, Martin et al. [14] leveraged
the spatiotemporal features of tweets to assess the responses on Hurricane
Matthew, and Yeo et al. [22] analyzed Twitter communications on 2016 Southern
Louisiana flood recovery process. Cheng et al. [2] studied the problem of
tracking disaster events through streaming data from Twitter. They proposed an
online matrix factorization-based method to identify topics across time
periods.

As people share their emotions [11] and opinions [3] related to the COVID-19
pandemic on social media (e.g. Twitter, Facebook), a number of studies have
also been conducted to explore the use of social media during COVID-19
epidemic. Although numerous researches have focused on online discourses [19,
6, 21], a few recent studies have tried to disclose the impact of various
offline data delivered by traditional news sources. Crupi et al. [4] showed
the change of topics regarding COVID-19 vaccination in Italian user
communities and the polarization of vaccine stances with common and distinct
topics. Another work by Poddar et al. [16] identified distinct topics in
online data and investigated how vaccine-related discourse has changed between
pre and during COVID times. A few studies have utilized geo-tagged data on
social media to map online data to real-world events and statistics [18, 5].
Others tried to uncover the association between vaccine hesitancy and online
user behavior [20, 13] or online misinformation [12, 15].

Our work is different than previous studies as we study the novel problem of
linking COVID-19-related offline and online data. The aim of this problem is
to examine how common and distinct topics around offline and online data
change over time, and how offline events affect online data.

##  3 Proposed Method

In this section, we first review the basic NMF model and its online variant
which is a modified matrix factorization for temporal data. Then, we discuss
the main components of the proposed method in detail. Inspired by Cheng et al.
[2], we employ the joint ONMF algorithm to extract the latent factors of
streaming data using matrix factorization. Considering a time span of 24
months, in every one-month time period, historical data is combined with newly
arriving data to identify common and distinct topics. We perform the joint
ONMF on COVID-19 tweets and its related news data to extract the topics.

###  3.1 Non-negative Matrix Factorization (NMF)

NMF [10] is a popular non-negative matrix decomposition algorithm which is
widely used for analysis of multivariate data such as latent representation of
text. Given a non-negative data matrix $\mathbf{V}$, NMF factorizes
$\mathbf{V}$ into two lower rank non-negative matrices $\mathbf{W}$ and
$\mathbf{H}$ as $\mathbf{V}\approx\mathbf{W}\mathbf{H}$. Matrix
$\mathbf{V}\in\mathbb{R}^{m\times n}$ indicates $m$ data, each having $n$
attributes. Matrix $\mathbf{H}$ is a feature matrix with each row represents a
latent factor. $\mathbf{W}\in\mathbb{R}^{m\times k}$ is a coefficient matrix,
reflecting the association weights between the data samples and the extracted
latent factors. NMF is solved by minimizing the distance between the original
matrix $\mathbf{V}$ and the reconstructed matrix
$\mathbf{V^{\prime}}=\mathbf{W}\mathbf{H}$.

The NMF method assumes both latent factors and input data are static and do
not change over time. Thus, we cannot directly apply this method for our
problem. One simple solution to apply NMF on streaming data is to feed NMF
with updated data matrix as new subset of the data joins at each time step
$t$. However, this approach underutilizes the previous NMF results in the
upcoming factorization results that is calculated for the newly collected data
[1]. Cao et al. [1] proposes Online-NMF (ONMF) to decompose matrix
$\mathbf{V}_{t+1}$ based on a new batch of data $\mathbf{U}$ at time step
$t+1$ and the previous matrix $\mathbf{V}_{t}$. Formally, the ONMF problem is
defined as $\mathbf{V}_{t+1}=\begin{pmatrix}\mathbf{V}_{t}\\\
\mathbf{U}\end{pmatrix}\approx\mathbf{W}_{t}\mathbf{H}_{t}$, which is solved
by minimizing the mean squared error between $\mathbf{V}$ and $\mathbf{WH}$
[1].

###  3.2 Linking Offline and Online Data

ONMF can efficiently update the latent factors (i.e., topics) in streaming
data. However, it is hard to see the global semantic changes in the generated
topics. To solve this problem, inspired by [2], we use a modified ONMF-based
approach, joint ONMF, to discover commonness and distinctiveness in topics.
Assuming there are $k$ topics in the given documents, joint ONMF identifies
$k_{c}$ and $k_{d}=k-k_{c}$ which are defined as the number of common topics
and distinct topics, respectively. To alleviate the massive computational cost
proportional to the growing data matrix $\mathbf{V}$, this method uses one of
the decomposed low rank matrices $\mathbf{H}_{t}$ to derive topics. Since
$\mathbf{H}_{t}$ is fixed at time step $t$, we perform the linear
transformation $\mathbf{H}^{*}\approx\mathbf{L}^{*}\mathbf{H}_{t}$ to
dynamically tune the dependency between $\mathbf{H}_{t}$ and newly arrived
data $\mathbf{U}$. In this transformation,
$(\mathbf{L}^{*}\in\mathbb{R}^{k\times k})$ is a low rank matrix that is used
to adjust the dependency between $\mathbf{H}_{t}$ and $\mathbf{U}$.

Joint ONMF aims to minimize the distance between common topics in
$\mathbf{H}^{*}$ and $\mathbf{H}_{U}$ and maximize the distance between
distinct topics in $\mathbf{H}^{*}$ and $\mathbf{H}_{U}$. Finally, we use the
following objective function to find common and distinctive topics:

$\displaystyle\min_{\mathbf{W}_{U},\mathbf{H}_{U},\mathbf{H}^{*},\mathbf{L}^{*}}\frac{1}{2}\|\mathbf{H}^{*}-\mathbf{L}^{*}\mathbf{H}_{t}\|^{2}_{F}+$
$\displaystyle\frac{1}{2}\|\mathbf{U}-\mathbf{W}_{U_{c}}\mathbf{H}_{U_{c}}-\mathbf{W}_{U_{d}}\mathbf{H}_{U_{d}}\|^{2}_{F}$
$\displaystyle+$ $\displaystyle\alpha
f_{c}(\mathbf{H}^{*}_{c},\mathbf{H}_{U_{c}})+\beta
f_{d}(\mathbf{H}^{*}_{d},\mathbf{H}_{U_{d}}),$ (1)

where the second term denotes the reconstruction error on matrix $\mathbf{U}$,
which is factorized into $\mathbf{W}_{U}$ and $\mathbf{H}_{U}$
$(\mathbf{U}\approx\mathbf{W}_{U}\mathbf{H}_{U})$. $f_{c}$ and $f_{d}$ are the
commonness score and distinctiveness score described in section 4. We leverage
parameters $\alpha$ and $\beta$ to find the optimal balance between three
different objectives: the accuracy of matrix reconstruction, common, and
distinct topics.

This method identifies the common and distinct topics over time within one
dataset. To examine the relationship between the datasets with different
sources, we apply the method on both datasets simultaneously. Specifically, we
split each dataset by months as delimiters and concatenate the subsets of
tweets and news pieces one after another along the timeline. (e.g., Tweets of
January 2020 $||$ News of January 2020 $||$ Tweets of February 2020 $||$ News
of February 2020 $\cdots$). Starting with the Tweets of January 2020, the
model outputs the common topics between each month (i.e., topics that appear
both in $t$ and $t+1$, considering $t$ as month) and distinct topics (i.e.,
two unique topics at time before $t$ and time $t+1$) whenever new data arrives
at every time stamp $t$. The result distinct topics could give an insight on
the main theme and statements raised within a group of people on Twitter,
along with the important events recorded by the press happened in the specific
time period. On the contrary, the output of common topics implies the tweets’
topics are closely related to the news pieces.

##  4 Experiments

We conduct qualitative and quantitative evaluation on the performance of the
joint ONMF model for extracting topics and gaining insights from both online
and offline data. We aim to answer two research questions: (Q1) How well joint
ONMF extracts common and distinct topics across online and offline data? (Q2)
How does the offline data impact the online data? To answer Q1, we compare the
performance of joint ONMF with baselines based on various evaluation metrics.
To answer Q2, we design experiments to measure the distinctiveness of online
and offline data.

###  4.1 Experimental Design

We use the CoVaxNet dataset111https://github.com/jiangbohan/CoVaxNet that
includes $4,790$ offline news articles and $1,831,220$ online tweets within a
time span of two years from January 2020 till December 2021 [7]. We further
use Newspaper3k222https://github.com/codelucas/newspaper to extract title,
body, and summary from news pieces. The online data is first filtered to
select $20K$ representative random tweets. Finally, we extract the initial
textual features using TF-IDF scores from both datasets. A large TF-IDF value
indicates that the term can better distinguishes documents from each other.

To design experiments, we start with online data from January 2020 followed by
offline data from the same time period. Assuming the topics of tweets vary in
accordance with newly uploaded news pieces within certain period, we assign
tweets and news by turns monthly. We choose $k_{c}=2$ and $k_{d}=3$. The
optimization parameters $\alpha$ and $\beta$ of subsection 3.2 are set to
$1,000$ and $0.1$, respectively. Finally, the step size of the joint ONMF
method is set to $100$.

To evaluate the effectiveness of joint ONMF, we consider the three baselines.
Standard NMF (SNMF) [9] approach repeatably calculates the latent factors
using the entire dataset whenever a new batch of data is received. ONMF [1]
method is a modified variant of NMF that focuses on handling large-scale and
streaming data. Finally, Pseudo-Deflation [8] approach uses a variant of ONMF,
to identify common and distinct topics between two documents. For fair
comparison between joint ONMF and the baselines, we evaluate the retrieved
topics within the same time period and hyperparameters.

(a) CScore

(b) DScore

(c) RE

(d) Training Time

Figure 2: Performance comparison of different methods

###  4.2 Experimental Results

Following Kim et al. [8], we evaluate the performance of joint ONMF using
commonness score, distinctiveness score, and reconstruction error:

Commonness Score (CScore) stands for the similarity between the $k$ common
topics at time $t$ and $t+1$ for every time stamp, calculated as
$\text{CScore}=\frac{1}{k_{c}}\|\mathbf{H}^{*}_{c}-\mathbf{H}_{\mathbf{U}_{c}}\|^{2}_{F}$.
A smaller CScore indicates better efficiency of the method on extracting
common topics. Difference Score (DScore) uses average symmetric Kullback-
Leibler (KL) divergence between all the distinct topic pairs. A large DScore
indicates the distribution of the obtained distinct topics are discrete.

$\displaystyle\text{DScore}=\frac{1}{2k^{2}_{d}}\sum^{k_{d}}_{i}\sum^{k_{d}}_{j}[$
$\displaystyle
h^{*i}_{d}\log(h^{i}_{d})^{T}+h^{i}_{\mathbf{U}_{d}}\log(h^{i}_{\mathbf{U}_{d}})^{T}-$
$\displaystyle-h^{*i}_{d}\log(h^{j}_{d})^{T}-h^{j}_{\mathbf{U}_{d}}\log(h^{*i}_{\mathbf{U}_{d}})^{T}]$
(2)

Reconstruction Error (RE) estimates the loss of the NMF on the newly arriving
data $\mathbf{U}$ at each time stamp $t$ based on
$\min_{\mathbf{W},\mathbf{H}}\|\mathbf{V}-\mathbf{W}\mathbf{H}\|$. Smaller RE
indicates the initial data matrix $\mathbf{U}$ can be reconstructed
accurately.

(a) CScore

(b) DScore

(c) RE

Figure 3: Performance comparison with different data metrics of offline data

We answer the two questions Q1 and Q2 using quantitative and qualitative
experiments to evaluate the quality of common and distinct topics in the
integrated dataset, and study the effect of offline data on the online data.

Quantitative Results:

Performance (Q1). To calculate CScore and DScore, we directly apply them on
joint ONMF and Pseudo-Deflation as they explicitly state which topics are the
common and distinct. For other baseline models we assign topics with the
smallest CScore as common topics $k_{c}$, while the others are treated as
distinct ones $k_{d}$. Figure 2 illustrates the evaluation results including
computation time of each method. From these results, we conclude that the
joint ONMF method outperforms ONMF and SNMF with a large gap in CScore and
DScore. Pseudo-Deflation has comparable CScore, however, the DScore is
significantly lower than joint ONMF. DScores of all four methods fluctuates
severely in later time periods due to the noisy and unbalanced nature of the
dataset. Among the baselines, Pseudo-Deflation presents remarkably smaller RE.
Besides Pseudo-Deflation, SNMF can often obtain the best RE because the goal
of SNMF focuses on minimizing RE while other online topic models seek to
balance between RE and computational efficiency. The joint ONMF method
presents competitive RE compared to SNMF and ONMF as it separately conducts
matrix factorization on data matrix $\mathbf{U}$.

Considering all three evaluation metrics and the computational cost, the joint
ONMF approach attains the most optimal balance among the baselines. For
instance, it can have the smallest CScore and the largest DScore, meanwhile
achieving negligible reconstruction error similar to SNMF. 2(d) shows that
joint ONMF is slightly slower than two baselines yet maintains consistent
training time, making it comparatively computationally efficient.

Impact (Q2). According to Figure 3, one interesting observation is that the
amount of offline information inconsistently influence the performance of
joint ONMF over each measurements. As mentioned earlier, we measure the
performance on {‘title’, ‘summary’, ‘body’} of each news piece, assuming that
‘title’ text has the smallest, and the ‘body’ text has the largest quantity of
information. Results show that given more information, the model can better
reconstruct original data but worse in distinguishing distinct topics.

Table 1: Examples of topics from Jan 1st to Mar 31st, 2020. $D_{T}$ and $D_{N}$ show distinct topics of tweet and news title, respectively. $C$ shows the common topics of between offline and online data. Jan 2020 | Feb 2020 | Mar 2020
---|---|---
$D_{T}$ |  $C$ |  $D_{N}$ |  $D_{T}$ |  $C$ |  $D_{N}$ |  $D_{T}$ |  $C$ |  $D_{N}$
misinfor- mation |  oil |  antigen |  deadly |  flu |  Remdesivir |  misinfor- mation |  patent |  immunity
children |  Japan |  China |  fund |  antivaccine movement |  influenza season |  worry |  mask wear |  Donald Trump
thank |  religious exemption |  HPV cancer |  happy |  pediatric |  South Korea |  child killer |  candidates |  CDC
kill |  Ebola |  Donald Trump |  thank |  conspiracy |  afford insurance |  propaganda |  INOVIO (Pharmaceuticals) |  adjuvants
amazing |  Johnson & Johnson |  tuber- culosis |  autism |  Pneumonia |  peptid |  prevent disease |  egg-based vaccine |  Summit

Qualitative Study:

Performance (Q1). We perform an extensive qualitative study on the validity of
the identified topics. Table 1 shows top five ranked topics discovered by the
joint ONMF method on the integrated online and offline data during the first
three time periods of the pandemic. The extracted common topics are closely
related to temporal and/or newly rising events (e.g., Ebola, Remdesivir, egg-
based vaccine) which are compelling and provocative as well (i.e., topics
discussed online change over time and selective real-world events). Meanwhile,
distinct topics consist of specific keywords relevant to COVID-19 or medical
events. Note that according to Table 1, exceptionally large portion of
emotional words account for the common topics of the online data.

Impact (Q2). We take a deeper analysis on the common topics to observe how the
two datasets are linked. For example, in Table 2, we sample tweets and news
headlines containing the topic keyword oil. It is worth mentioning that the
users not only talk about the facts and opinions, but also they expose their
strong emotions (e.g., i will laugh, lol) and groundless claims (e.g., oil
companies just hid their own research data) that could falsely lead to
misinformation. We could find similar characteristics in multiple other tweet
samples. This implies that the offline data (i.e., local events and news) play
a substantial role in the shift of online behavior and potential biases of
users.

Table 2: Examples of different sources from common topic  Topic | Source |  Examples
---|---|---
oil | News Articles |  \- Russia halts oil to Belarus, but transit to Europe still flowing (Jan 3, 2020)
Tweets |  \- Lol oh yes definitely max single handedly created the need smh he is just one of the … antivaxxers of the oil world (Jan 4, 2020)
\- The oil companies just hid their own research data now all the sceptics
have left like flat earthers … (Jan 5, 2020)
\- I will laugh at anyone who sneers at China for this given that America is a
rapidly decaying shithole filled with … climate denialism antivaxxers and oil
companies that own multiple politicians (Jan 31, 2020)
trial | News Articles |  \- Coronavirus Vaccine: NIH Goes Straight to Human Trials (Mar 16, 2020)
Tweets |  \- Let’s go science NIH clinical trial of investigational vaccine for covid 19 begins … (Mar 16, 2020)
\- Why i am volunteering to get the fastest startup of a vaccine trial ever …
(Mar 31, 2020)
delta | News Articles |  \- Pfizer says COVID vaccine is highly effective against Delta variant (Jun 24, 2021)
\- Moderna Says Studies Show Its Vaccine Is Effective Against The Delta
Variant (Jun 30, 2021)
Tweets |  \- This is what I’ve been saying the delta variant is growing and now there is a new lamda variant while we can’t live in fear over every variant that pops up we will still need to be cautious until we learn what they do and whom they affect (July 8, 2021)
\- Get vaccinated or get sick but do not punish those who are vaccinated I am
fed up with restrictions on us because of and citing rise of delta variant Los
Angeles reports 165 percent increase in covid cases (July 9, 2021)

##  5 Conclusion and Future Work

In this paper, we study the novel problem of integrating online and offline
COVID-19-related data to seek common and distinct topics across them. Finding
the commonness and distinctiveness between online and offline data is an
important research problem because it provides an effective way for
researchers to study how users’ opinions change over time and how offline news
sources affect online social media discussions. Furthermore, studying this
problem paves the way for internet users to gather information on
COVID-19-related topics and avoid possible misinformation effectively. To this
end, we use the joint ONMF method to compute common and distinct topics
efficiently. We design qualitative and quantitative experiments to measure the
performance of joint ONMF and provide valuable insights into how offline data
affects online data. For future work, we will study methods that can automate
the process of finding interesting points of view and user attributes
contributed to online discussions about offline events.

## Acknowledgement

This work was supported by the Office of Naval Research under Award No.
N00014-21-1-4002. Opinions, interpretations, conclusions, and recommendations
are those of the authors.

## References

  * [1] Cao, B., et al.: Detect and track latent factors with online nonnegative matrix factorization. In: IJCAI (2007)
  * [2] Cheng, L., et al.: Tracking disaster footprints with social streaming data. In: IAAI (2020)
  * [3] Cotfas, L.A., et al.: The longest month: analyzing covid-19 vaccination opinions dynamics from tweets in the month following the first vaccine announcement. IEEE Access (2021)
  * [4] Crupi, G., et al.: Echoes through time: Evolution of the italian covid-19 vaccination debate. In: ICWSM (2022)
  * [5] Feng, S., Kirkley, A.: Integrating online and offline data for crisis management: Online geolocalized emotion, policy response, and local mobility during the covid crisis. Scientific Reports (2021)
  * [6] Glandt, K., et al.: Stance detection in covid-19 tweets. In: IJCNLP (2021)
  * [7] Jiang, B., Sheth, P., Li, B., Liu, H.: Covaxnet: An online-offline data repository for covid-19 vaccine hesitancy research. arXiv preprint arXiv:2207.01505 (2022)
  * [8] Kim, H., et al.: Simultaneous discovery of common and discriminative topics via joint nonnegative matrix factorization. In: KDD (2015)
  * [9] Lee, D., Seung, H.S.: Algorithms for non-negative matrix factorization. NeurIPS (2000)
  * [10] Lee, D.D., et al.: Learning the parts of objects by non-negative matrix factorization. Nature (1999)
  * [11] Lwin, M.O., et al.: Global sentiments surrounding the covid-19 pandemic on twitter: analysis of twitter trends. JMIR public health and surveillance (2020)
  * [12] Lyu, H., et al.: Misinformation versus facts: Understanding the influence of news regarding covid-19 vaccines on vaccine uptake. Health Data Science (2022)
  * [13] MacDonald, N.E., et al.: Vaccine hesitancy: Definition, scope and determinants. Vaccine (2015)
  * [14] Martín, Y., et al.: Leveraging twitter to gauge evacuation compliance: Spatiotemporal analysis of hurricane matthew. PLoS one (2017)
  * [15] Pierri, F., et al.: Online misinformation is linked to early covid-19 vaccination hesitancy and refusal. Scientific reports (2022)
  * [16] Poddar, S., et al.: Winds of change: Impact of covid-19 on vaccine-related opinions of twitter users. In: ICWSM (2022)
  * [17] Reynolds, B., et al.: Crisis and emergency risk communication (2002)
  * [18] Shen, C., et al.: Using reports of symptoms and diagnoses on social media to predict covid-19 case counts in mainland china: Observational infoveillance study. JMIR (2020)
  * [19] Wicke, P., Bolognesi, M.M.: Covid-19 discourse on twitter: How the topics, sentiments, subjectivity, and figurative frames changed over time. Frontiers in Communication (2021)
  * [20] Wilson, S.L., Wiysonge, C.: Social media and vaccine hesitancy. BMJ global health (2020)
  * [21] Xiong, Z., et al.: Social media opinions on working from home in the united states during the covid-19 pandemic: Observational study. JMIR medical informatics (2021)
  * [22] Yeo, J., et al.: Disaster recovery communication in the digital era: Social media and the 2016 southern louisiana flood. Risk analysis (2020)


# A Simple Characterization of Supply Correspondences††thanks:
_Krishnamoorthy:_ Carnegie Mellon University, Pittsburgh, USA;
<EMAIL_ADDRESS>_Kushnir:_ Carnegie Mellon University, Tepper School
of Business, Pittsburgh, USA<EMAIL_ADDRESS>We are very thankful to
Paul Milgrom for sharing lecture notes on producer choice theory by
LevinMilgromSegal2016. They inspired this project.

Vinod Krishnamoorthy  Alexey Kushnir

(May 2022)

###### Abstract

We prove that supply correspondences are characterized by two properties: the
law of supply and being homogeneous of degree zero.

JEL classification: D21, D24
Keywords: supply correspondence, supply function, rationalizability, the law
of supply, monotonicity, cyclic-monotonicity, homogeneous of degree zero.

##  1 Introduction

We prove that supply correspondences are characterized by two simple
properties: the law of supply and being homogeneous of degree zero. These two
properties are both basic properties of profit maximization behavior: The
first one captures the intuition that the supply decisions should change in
the same direction as prices change. The second property states that if all
prices of inputs and outputs are proportionally changed, then the firm supply
decision should not change (see Section 5.C in mas1995microeconomic).

To address the question of when a given supply correspondence is consistent
with profit maximization behavior, we consider two concepts of
rationalizability. A supply correspondence is rationalizable if there exists a
convex closed production set such that for every price vector each supply
decision is profit maximizing. A supply correspondence is strongly
rationalizable if it consists of all possible maximizers. Our first result
shows that if a supply correspondence satisfies the law of supply and is
homogeneous of degree zero, then it is rationalizable. To obtain the strong
form of rationalizability, we strengthen the law of supply property. The law
of supply is equivalent to a supply correspondence being monotone. Consider
the set of all monotone correspondences together with their graphs. A maximal
monotone correspondence is a monotone correspondence with the graph that
cannot be a strict subset of the graph of any other monotone correspondence;
that is, its graph is maximal by set inclusion. We show that for a supply
correspondence to be strongly rationalizable, it is necessary and sufficient
for it to be maximal monotone and homogeneous of degree zero.

The previous literature has provided several characterizations of
rationalizable supply functions.111See samuelson1948foundations,
hanoch1972testing, Varian1984 for earlier important contributions to
production analysis. Most of these characterizations make strong assumptions
about the differentiability of supply function. For example, Proposition 8 in
LevinMilgromSegal2016 states that a continuously differentiable supply
function $y$ is rationalizable if and only if it satisfies $Dy(p)p=0$ and its
Jacobian $Dy(p)$ is symmetric, positive semidefinite.222For a similar
statement, see Proposition 7.9 in kreps2013microeconomic. jehle2011advanced in
Chapter 3 reference this result as integrability theorem supply functions,
which is parallel to a similar result in consumer choice theory. See also
mas1995microeconomic and varian92. The condition $Dy(p)p=0$ is nothing but
Euler’s law applied to a homogeneous of degree zero supply function. The
condition on Jacobian is an analog of the cyclic-monotonicity condition. Our
characterization applies to supply correspondences. When supply correspondence
satisfies the law of supply, we show that the cyclic-monotonicity condition
could be dropped or reduced to the maximal monotonicity condition, depending
on what rationalizability concept one is interested in.333An important result
in convex analysis is that if a correspondence is maximal monotone and
cyclically monotone, then it is the subdifferential of some proper convex
lower semicontinuous function (i.e., the profit function) (see
rockafellar2015convex; phelps1997lectures). In relation to this result, we
show that the cyclic-monotonicity condition can be dropped if the function is
homogeneous of degree zero. Homogeneous of degree zero is a strong condition,
but it is very natural in the context of production theory. Both conditions
are easy to verify. Also, our characterizations does not require any
assumptions of differentiablity on the supply function.

Our characterization of supply correspondences is also closely related to the
characterization of combinatorial demand correspondences in the recent paper
by chambers2018characterization. However, we note several important
differences. First, our characterization does not require the image of supply
correspondence to be finite (combinatorial). Second, we have different
properties: we characterize supply correspondences in terms of the law of
supply and homogeneity of degree zero, whereas chambers2018characterization
characterize combinatorial demand correspondences in terms of the law of
demand and upper hemicontinuity. While both the law of demand and the law of
supply are two versions of the monotonicity condition with opposite signs, the
homogeneity of degree zero condition is absent in consumer choice theory, as
the maximization objective is not linear. We also illustrate in Section
LABEL:sec:discussion why the maximal monotonicity condition cannot be replaced
with the upper hemicontinuity condition when the image of supply
correspondence is not finite.

##  2 Notation

Let us consider an economy with $N\geq 2$ commodities. A typical production
plan is a vector $z=(z_{1},...,z_{N})\in\mathbb{R}^{N}$, where an output has
$z_{n}>0$ and an input has $z_{n}<0$. The vector of prices is denoted as
$p\in\mathbb{R}^{N}$.444We could alternatively assume that the set of possible
prices is any cone containing the origin. In particular, we can assume that
prices cannot be negative; that is, $p\in\mathbb{R}^{N}_{+}$ (see the
discussion in Section LABEL:sec:discussion). We say that for
$p,p^{\prime}\in\mathbb{R}^{N}$, $p\leq p^{\prime}$ if each coordinate of $p$
is smaller than the corresponding coordinate of $p^{\prime}$.

A supply correspondence is a mapping
$y:\mathbb{R}^{N}\rightrightarrows\mathbb{R}^{N}$, giving for each price
vector $p$ a set of possible production plans. If $y$ is single-valued, we
refer to it as a supply function. Let
$\textrm{Im}(y)=\\{z\in\mathbb{R}^{n}\,|\,\exists p\in\mathbb{R}^{N}\text{
s.t. }z\in y(p)\\}$. We consider two concepts of rationalizability.

######  Definition 1.

Supply correspondence $y:\mathbb{R}^{N}\rightrightarrows\mathbb{R}^{N}$ is
rationalizable if there exists a convex and closed production set $Y$ such
that for all $p\in\mathbb{R}^{N}$,
$y(p)\subseteq\\{z\in\mathbb{R}^{N}\,|\,p\cdot z=\sup_{z^{\prime}\in Y}p\cdot
z^{\prime}\\}.$

We note that the requirement for production set $Y$ to be convex and closed is
innocuous, because if we find some production set $Y$ to satisfy the
maximization condition, the closure of its convex hull will rationalize supply
correspondence $y$ as well.