hypothesisHypothesis claimClaim Double precision is not necessary for LSQR H.
Li, G. Tan and T. Zhao

# Double precision is not necessary for LSQR for solving discrete linear ill-
posed problems ††thanks: This work of the second author was supported in part
by the National Natural Science Foundation of China under Grant No. 62032023.

Haibo Li Corresponding author. Computing System Optimization Lab, Huawei
Technologies, Beijing 100094, China, and Institute of Computing Technology,
Chinese Academy of Sciences, Beijing 100190, China<EMAIL_ADDRESS>Guangming Tan Institute of Computing Technology, Chinese Academy of Sciences,
Beijing 100190, China<EMAIL_ADDRESS><EMAIL_ADDRESS>Tong Zhao
33footnotemark: 3

###### Abstract

The growing availability and usage of low precision floating point formats
attracts many interests of developing lower or mixed precision algorithms for
scientific computing problems. In this paper we investigate the possibility of
exploiting lower precision computing in LSQR for solving discrete linear ill-
posed problems. We analyze the choice of proper computing precisions in the
two main parts of LSQR, including the construction of Lanczos vectors and
updating procedure of iterative solutions. We show that, under some mild
conditions, the Lanczos vectors can be computed using single precision without
loss of any accuracy of the final regularized solution as long as the noise
level is not extremely small. We also show that the most time consuming part
for updating iterative solutions can be performed using single precision
without sacrificing any accuracy. The results indicate that the most time
consuming parts of the algorithm can be implemented using single precision,
and thus the performance of LSQR for solving discrete linear ill-posed
problems can be significantly enhanced. Numerical experiments are made for
testing the single precision variants of LSQR and confirming our results.

###### keywords:

lower precision, linear ill-posed problem, regularization, LSQR, roundoff
unit, semi-convergence

65F10, 65F22, 65F50, 65G50

##  1 Introduction

Although for most traditional scientific computing problems, computations are
carried out with double precision (64-bit) rather than lower precisions such
as single (32-bit) or half (16-bit) precision, on modern computing
architectures, the performance of 32-bit operations is often at least twice as
fast as that of 64-bit operations, which stimulates the trial of using lower
precision floating point formats in an algorithm. To exploit this computation
power without sacrificing accuracy of the final result, numerical algorithms
have to be designed that use lower or mixed precision formats. By using a
combination of 64-bit and 32-bit (even 16-bit) floating point arithmetic, the
performance of many numerical algorithms can be significantly enhanced while
maintaining the 64-bit accuracy of the resulting solution. New mixed precision
variants of many numerical linear algebra algorithms have been recently
proposed, such as matrix multiplications [1, 5], LU and QR matrix
factorizations [28, 39], Krylov solvers [15, 7] and many others [22, 2].

In this paper, we investigate how to exploit lower precision computing for
solving discrete linear ill-posed problems. This type of problems typically
arise in connection with the numerical solution of inverse problems that
appear in various applications of science and engineering, such as image
deblurring, geophysics, computerized tomography and many others; see e.g.,
[31, 24, 20, 35]. A basic linear inverse problem leads to a discrete linear
system of the form

(1) $Ax=b,\ \ \ \ b=Ax_{ex}+e,$

where the matrix $A\in\mathbb{R}^{m\times n}$ with $m\geq n$ without loss of
generality, and the right-hand side $b$ is a perturbed version of the unknown
exact observation $b_{ex}=Ax_{ex}$. In this paper we suppose $e$ is a Gaussian
white noise. The problem is ill-posed in the sense that $A$ is extremely ill-
conditioned with its singular values decaying gradually towards zero without
any noticeable gap, which leads to that the naive solution
$x_{nai}=A^{{\dagger}}b$ of Eq. 1 is a poor approximation to the exact
solution $x_{ex}=A^{{\dagger}}b_{ex}$, where $``{\dagger}"$ denotes the Moore-
Penrose inverse of a matrix. Therefore, some forms of regularization must be
used to deal with the noise $e$ in order to extract a good approximation to
$x_{ex}$.

One of the popular regularization techniques is the Tikhonov regularization
[36], in which a quadratic penalty is added to the objective function:

$x_{\lambda}=\arg\min_{x\in\mathbb{R}^{n}}\\{\lVert
Ax-b\lVert^{2}+\lambda\lVert Lx\lVert^{2}\\},$

where $\lambda>0$ is the regularization parameter and $L\in\mathbb{R}^{p\times
n}$ is the regularization matrix. Throughout the rest of the paper $\|\cdot\|$
always denotes either the vector or matrix 2-norm. The proper choice of $L$
depends on the particular application, which should be chosen to yield a
regularized solution with some known desired features of $x_{ex}$. A suitable
value of $\lambda$ should have a good balance between the fidelity term and
the regularization term, only in which case we can get a regularized solution
that is a good approximation to $x_{ex}$. Although many types of
regularization parameter choice rules have been proposed, such as discrepancy
principle (DP) [30], unbiased predictive risk estimator [37], generalized
cross validation [14] and L-curve method [16], it is often computationally
very expensive to choose a suitable $\lambda$ for large scale problems, since
many different values of $\lambda$ must be tried to get $x_{\lambda}$.

For large scale ill-posed problems, iterative regularization method is the
soundest choice. For standard-form regularization with $L=I_{n}$, Krylov
subspace based methods [3, 12] are the most commonly used. The method projects
(1) onto a sequence of lower-dimensional Krylov subspaces, and then solves the
projected small scale problems, where the iteration number plays the role of
regularization parameter [17, 19]. This approach usually exhibits semi-
convergence: at early iterations solutions converge to $x_{ex}$ while
afterwards the noise $e$ starts to deteriorate the solutions so that they
start to diverge from $x_{ex}$ and instead converge to $x_{nai}$. In order to
get a good regularized solution, the iteration must be stopped early properly
by using a regularization parameter choice rule [6, 33]. The semi-convergence
behavior can be mitigated by using a hybrid method, which applies a standard
regularization technique, such as Tikhonov regularization or truncated SVD, to
the projected problem at each iteration [25, 9, 34].

In this paper, we focus on the LSQR algorithm for iteratively solving large
scale discrete linear ill-posed problems, which is a based on the Lanczos
bidiagonalization [32]. The motivation for this work is to answer whether
lower precisions can be used in some parts of the algorithm while maintaining
the 64-bit accuracy of the regularized solutions. A crucial problem for LSQR
in finite precision arithmetic is the loss of orthogonality of Lanczos vectors
constructed by the Lanczos bidiagonalizatoin, which leads to slower and
irregular convergence of iterative solutions [26]. When using lower precision
storage and arithmetic, errors attributed to rounding errors may have a larger
influence on the total error, and the accuracy of the final result may be more
sensitive to rounding errors. Thus the LSQR is usually implemented with full
reorthogonalization of Lanczos vectors for solving ill-posed problems [23]. We
will investigate how to implement lower precision variants of LSQR with full
reorthogonalization to solve discrete linear ill-posed problems.

Since $b$ is contaminated by the noise $e$, we can never get a regularized
solution with error as small as $\|e\|$ [10, 17], and this error is usually
much bigger than the roundoff unit of double precision that is
$\mathbf{u}=2^{-53}$. This fact inspires us that double precision may be not
necessary for LSQR to compute a regularizes solution with the same accuracy as
the best regularized one. In this paper, we study lower precision computing
for the two main parts of LSQR, including the construction of Lanczos vectors
and updating procedure of iterative solutions. Under an ideal model describing
the ill-posed problem (1), our result can determine which level of computing
precision should be used for constructing Lanczos vectors with full
reorthogonalization, and it indicates that for not extremely small noise
levels we can exploit single precision for this part without loss of accuracy
of the final regularized solution. For the updating procedure part, we
theoretically show that, under a condition which can be almost always
satisfied, the updated solutions can be computed using single precision
without sacrificing any accuracy. The results show that the most time
consuming parts of LSQR for solving discrete linear ill-posed problems can be
implemented using single precision, and thus the performance can be
significantly enhanced. The results can guide us towards an implementation
that is both robust and efficient, and should be considered by application
developers for their practical problems.

The paper is organized as follows. We start with a brief review of LSQR for
discrete linear ill-posed problems and the IEEE 754 floating point standard in
Section 2. Next, in Section 3, we make an analysis to determine which level of
computing precision should be used for constructing Lanczos vectors with full
reorthogonalization. In Section 4, we analyze the possibility of using single
precision for updating iterative solutions, and also discuss the estimation of
semi-convergence point of LSQR in lower precision arithmetic. Some numerical
experimental results are presented in Section 5 and we conclude the paper in
Section 6.

##  2 Preliminaries

In this section, we review some basic knowledge of discrete linear ill-posed
problems and the LSQR regularization algorithm, we also review the IEEE 754
Standard floating point number system.

Suppose the singular value decomposition (SVD) of $A$ is

(2) $A=\hat{U}\left(\begin{array}[]{c}\Sigma\\\
\mathbf{0}\end{array}\right)\hat{V}^{T},$

where $\hat{U}=(\hat{u}_{1},\ldots,\hat{u}_{m})\in\mathbb{R}^{m\times m}$ and
$\hat{V}=(\hat{v}_{1},\ldots,\hat{v}_{n})\in\mathbb{R}^{n\times n}$ are
orthogonal, $\Sigma={\rm
diag}(\sigma_{1},\ldots,\sigma_{n})\in\mathbb{R}^{n\times n}$ with singular
values $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{n}>0$. The naive
solution to (1) is

$x_{nai}=\sum\limits_{i=1}^{n}\frac{u_{i}^{T}b}{\sigma_{i}}v_{i}=\sum\limits_{i=1}^{n}\frac{u_{i}^{T}b_{ex}}{\sigma_{i}}v_{i}+\sum\limits_{i=1}^{n}\frac{u_{i}^{T}e}{\sigma_{i}}v_{i}=x_{ex}+\sum\limits_{i=1}^{n}\frac{u_{i}^{T}e}{\sigma_{i}}v_{i}.$

Note that the second term in $x_{nai}$ is extremely large since $\sigma_{i}$
decay to zero, making $x_{nai}$ a meaningless solution. For discrete linear
ill-posed problems, the exact discrete Picard condition(DPC) plays a central
role, which says that the Fourier coefficients $|\hat{u}_{i}^{T}b_{ex}|$ on
average decay to zero faster than the corresponding $\sigma_{i}$ [17, 19]. It
is well known that any regularization is based on an underlying requirement
that the DPC for a given problem is satisfied, only under which can one
compute a useful regularized solution with some accuracy [17, 19]. Throughout
the rest of the paper, we assume that the DPC is satisfied and
$\|e\|<\|b_{ex}\|$ (otherwise all informations from $b_{ex}$ is lost in $b$).

A direct regularized method is the truncated SVD (TSVD) method, which forms
$x_{k}^{tsvd}$ by truncating the first $k$ components of $x_{nai}$
corresponding to large singular values:
$x_{k}^{tsvd}=\sum\limits_{i=1}^{k}\frac{u_{i}^{T}b}{\sigma_{i}}v_{i}$. Under
the assumption that $e$ is Gaussian, it is known from [17, 10] that there
exist $k^{*}$ such that $x_{k^{*}}^{tsvd}$ is a 2-norm filtering best possible
regularized solution of (1) when only deterministic 2-norm filtering
regularization methods are taken into account. For large scale ill-posed
problems, the LSQR algorithm with early stopping is the most common used
iterative regularization method, which is based on Lanczos bidiagonalization
of $A$ with starting vector $b$, as described in Algorithm 1.

Algorithm 1 $k$-step Lanczos bidiagonalization

1: Let $\beta_{1}=\|b\|,\ u_{1}=b/\beta_{1}$

2: $\alpha_{1}=\|A^{T}u_{1}\|,\ v_{1}=A^{T}u_{1}/\alpha_{1}$

3: for $j=1,2,\dots,k$ do

4: $p_{j}=Av_{j}-\alpha_{j}u_{j}$

5: $\beta_{j+1}=\|p_{j}\|$, $u_{j+1}=p_{j}/\beta_{j+1}$

6: $r_{j}=A^{T}u_{j+1}-\beta_{j+1}v_{j}$

7: $\alpha_{j+1}=\|r_{j}\|$, $v_{j+1}=r_{j}/\alpha_{j+1}$

8: end for

In exact arithmetic, the $k$-step Lanczos bidiagonalization produces a lower
bidiagonal matrix

$B_{k}=\begin{pmatrix}\alpha_{1}&&&\\\ \beta_{2}&\alpha_{2}&&\\\
&\beta_{3}&\ddots&\\\ &&\ddots&\alpha_{k}\\\
&&&\beta_{k+1}\end{pmatrix}\in\mathbb{R}^{(k+1)\times k},$

and two groups of Lanczos vectors $\\{u_{1},\cdots,u_{k+1}\\}$ and
$\\{v_{1},\cdots,v_{k+1}\\}$ that are orthonormal bases of Krylov subspaces
$\mathcal{K}_{k+1}(AA^{T},b)$ and $\mathcal{K}_{k+1}(A^{T}A,A^{T}b)$,
respectively. The $k$-step Lanczos bidiagonalization can be written in the
matrix form

(3) $\displaystyle U_{k+1}(\beta_{1}e_{1}^{(k+1)})=b,$ (4) $\displaystyle
AV_{k}=U_{k+1}B_{k},$ (5) $\displaystyle
A^{T}U_{k+1}=V_{k}B_{k}^{T}+\alpha_{k+1}v_{k+1}(e_{k+1}^{(k+1)})^{T},$

where $e_{i}^{(l)}$ denotes the $i$-th canonical basis vector of
$\mathbb{R}^{l}$, and $U_{k+1}=(u_{1},\dots,u_{k+1})$ and
$V_{k}=(v_{1},\dots,v_{k})$ are two orthonormal matrices. The LSQR for (1) is
mathematically equivalent to the conjugate gradient (CG) method applied to the
normal equation of $\min_{x\in\mathbb{R}^{n}}\|Ax-b\|$, i.e. $A^{T}Ax=A^{T}b$,
which seeks approximations to $x_{ex}$ from the $k$ dimensional Krylov
subspace $\mathcal{K}_{k}(A^{T}A,A^{T}b)=\mathcal{R}(V_{k})$ starting with
$k=1$ onwards, and the iteration should be terminated at a proper step to
avoid semi-convergence to get the best possible regularized solution [3]. At
the $k$-th step, by (3) and (4) we have

$\displaystyle\ \ \ \ \min\limits_{x=V_{k}y}\lVert Ax-b\lVert$
$\displaystyle=\min\limits_{y\in\mathbb{R}^{k}}\lVert
AV_{k}y-\beta_{1}U_{k+1}e_{1}^{(k+1)}\lVert$
$\displaystyle=\min\limits_{y\in\mathbb{R}^{k}}\lVert
U_{k+1}B_{k}y-\beta_{1}U_{k+1}e_{1}^{(k+1)}\lVert$
$\displaystyle=\min\limits_{y\in\mathbb{R}^{k}}\lVert
B_{k}y-\beta_{1}e_{1}^{(k+1)}\lVert\ .$

and thus the $k$-step LSQR solution is

(6) $x_{k}=V_{k}y_{k},\ \
y_{k}=\mathrm{arg}\min\limits_{y\in\mathbb{R}^{k}}\lVert
B_{k}y-\beta_{1}e_{1}^{(k+1)}\lVert=B_{k}^{{\dagger}}(\beta_{1}e_{1}^{(k+1)})\
.$

We note that if $\alpha_{k+1}$ or $\beta_{k+1}$ is zero, then the iteration
terminates and $x_{k}=x_{nai}$ [32]. This case rarely happens in real
computations and we assume that the iteration does not terminate.

From the above description, the computation of $x_{k}$ can be divided into two
parts. The first part is the Lanczos bidiagonalization that generates two
orthonormal bases of Krylov subspaces $\mathcal{K}_{k+1}(AA^{T},b)$ and
$\mathcal{K}_{k+1}(A^{T}A,A^{T}b)$, respectively, while the second part is
solving the projected problem (6) to obtain $x_{k}$. In fact, in the practical
implementation, there is a recursive formula to update $x_{k+1}$ from $x_{k}$
without solving the projected least squares problems at each iteration. We
will analyze the choice of proper computing precisions for these two parts in
the following two sections.

In practical computational tasks, the accuracy of a computed result and the
time needed to complete the algorithm both heavily depend on the floating
point format used for storage and arithmetic operations. Here we review the
IEEE 754 Standard floating point number format, which is composed of a sign
bit, an exponent $\eta$, and a significand $t$:

$x=\pm\mu\times 2^{\eta-t},$

where $\mu$ is any integer in $[0,2^{t}-1]$ and $\eta$ is an integer in
$[\eta_{\text{min}},\eta_{\text{max}}]$. Roughly speaking, the length of the
exponent determines the value range of a floating point format, and the length
of the significand determines the relative accuracy of the format in that
range. A short analysis of floating point operations [21, Theorem 2.2] shows
that the relative error is controlled by the roundoff unit
$\mathbf{u}:=\frac{1}{2}\cdot 2^{1-t}$. Table 2.1 shows main parameters of the
three different floating point number formats.

Table 2.1: Parameters for various floating-point formats. “Range” denotes the order of magnitude of the smallest positive (subnormal) $x_{min}^{s}$ and smallest and largest positive normalized floating-point numbers. Type | Size | Range | Roundoff unit
---|---|---|---
| (bits) | $x_{min}^{s}$ | $x_{min}$ | $x_{max}$ | $\mathbf{u}$
half precision | 16 | $5.96\times 10^{-8}$ | $6.10\times 10^{-5}$ | $6.55\times 10^{4}$ | $4.88\times 10^{-4}$
single precision | 32 | $1.40\times 10^{-45}$ | $1.18\times 10^{-38}$ | $3.40\times 10^{38}$ | $5.96\times 10^{-8}$
double precision | 64 | $4.94\times 10^{-324}$ | $2.22\times 10^{-308}$ | $1.80\times 10^{308}$ | $1.11\times 10^{-16}$

In finite precision arithmetic, the Lanczos vectors $u_{i}$ and $v_{i}$
computed by the Lanczos bidiagonalization gradually lose their orthogonality,
and this may cause the LSQR to require significantly more iterations to
converge and make the propagation of noise during the iterations rather
irregular [29, 23]. Since it usually takes not too many iterations of the LSQR
to reach the semi-convergence point, it is usually implemented with full
reorthogonalization 111In full reorthogonalization, $u_{k}$ and $v_{k}$ are
reorthogonalized against all previous vectors $\\{u_{1},\dots,u_{k-1}\\}$ and
$\\{v_{1},\dots,v_{k-1}\\}$ as soon as they have been computed. This adds an
arithmetic cost of about $4(m+n)k^{2}$ flops, which is affordable if
$k\ll\min\\{m,n\\}$. to maintain stability of convergence. In the rest of the
paper, we investigate the LSQR implemented with full reorthogonalization in
finite precision arithmetic.

Define the orthogonality level of Lanczos vectors $\\{u_{1},\dots,u_{k}\\}$
and $\\{v_{1},\dots,v_{k}\\}$ as

$\mu_{k}=\lVert\mathbf{SUT}(I_{k}-U_{k}^{T}U_{k})\lVert,\ \ \ \
\nu_{k}=\lVert\mathbf{SUT}(I_{k}-V_{k}^{T}V_{k})\lVert,$

where $\mathbf{SUT(\cdot)}$ denotes the strictly upper triangular part of a
matrix. The following result has been established for the $k$-step Lanczos
bidiagonalization with reorthogonalization (not necessarily full
reorthogonalization) [27].

######  Theorem 2.1.

For the $k$-step Lanczos bidiagonalization with reorthogonalization, if
$\nu_{k+1}<1/2$ and $\mu_{k+1}<1/2$, then there exist two orthornormal
matrices
$\bar{U}_{k+1}=(\bar{u}_{1},\dots,\bar{u}_{k+1})\in\mathbb{R}^{m\times(k+1)}$
and $\bar{V}_{k+1}=(\bar{v}_{1},\dots,\bar{v}_{k+1})\in\mathbb{R}^{n\times k}$
such that

(7) $\displaystyle\bar{U}_{k+1}(\beta_{1}e_{1}^{(k+1)})=b+\delta_{b},$ (8)
$\displaystyle(A+E)\bar{V}_{k}=\bar{U}_{k+1}B_{k},$ (9)
$\displaystyle(A+E)^{T}\bar{U}_{k+1}=\bar{V}_{k}B_{k}^{T}+\alpha_{k+1}\bar{v}_{k+1}(e_{k+1}^{(k+1)})^{T},$

where $E$ and $\delta_{b}$ are perturbation matrix and vector, respectively.
We have error bounds

$\lVert\bar{U}_{k+1}-U_{k+1}\lVert\leq 2\mu_{k+1}+O(\mu_{k+1}^{2}),\ \ \
\lVert\bar{V}_{k+1}-V_{k+1}\lVert\leq\nu_{k+1}+O(\nu_{k+1}^{2}),$

and

$\lVert E\lVert=O(c(n,k)\lVert A\lVert(\mathbf{u}+\nu_{k+1}+\mu_{k+1})),\ \ \
\lVert\delta_{b}\lVert=O(\lVert b\lVert\mathbf{u}).$

where $c(n,k)$ is a moderately growing constant depends on $n$ and $k$.

By Theorem 2.1, for the Lanczos bidiagonalization with full
reorthogonalization, the orthogonality levels of $u_{i}$ and $v_{i}$ are kept
around $O(\mathbf{u})$, and thus we have

(10) $\displaystyle\lVert\bar{U}_{k+1}-U_{k+1}\|=O(\mathbf{u}),\ \ \
\lVert\bar{V}_{k+1}-V_{k+1}\lVert=O(\mathbf{u}),$ (11) $\displaystyle\lVert
E\lVert=O(c(n,k)\lVert A\lVert\mathbf{u}),\ \ \
\lVert\delta_{b}\lVert=O(\lVert b\lVert\mathbf{u}).$

This result will be used in the next section to determine a proper value of
roundoff unit $\mathbf{u}$ of the exploited computing precision, with which
the LSQR can compute a final regularized solution with the same accuracy as
the best regularized one.

##  3 Choice of computing precision for constructing Lanczos vectors

In this section, we investigate which level of computing precision should be
used for constructing Lanczos vectors with full reorthogonalization. Here we
assume that (6) is solved in exact arithmetic. We need to give an upper bound
on roundoff unit $\mathbf{u}$ such that for a computing precision with this
$\mathbf{u}$ the computed $x_{k}$ can achieve the same accuracy as the best
regularized LSQR solution to (1) obtained in exact arithmetic.

######  Theorem 3.1.

Follow the notations in Theorem 2.1. Suppose that the LSQR is implemented with
full reorthogonalization. Then the computed $x_{k}$ satisfies

(12) $\dfrac{\|x_{k}-\bar{x}_{k}\|}{\|\bar{x}_{k}\|}=O(\mathbf{u}).$

where $\bar{x}_{k}$ is the exact $k$-th LSQR solution to the perturbed problem

(13) $\min\limits_{x\in\mathbb{R}^{n}}\lVert(A+E)x-(b+\delta_{b})\lVert.$

######  Proof 3.2.

Since (6) is solved exactly, by (7) and (8) we have

$\displaystyle\begin{split}y_{k}&=\mathrm{arg}\min\limits_{y\in\mathbb{R}^{k}}\lVert
B_{k}y-\beta_{1}e_{1}^{(k+1)}\lVert=\mathrm{arg}\min\limits_{y\in\mathbb{R}^{k}}\lVert\bar{U}_{k+1}B_{k}y-\bar{U}_{k+1}\beta_{1}e_{1}^{(k+1)}\lVert\\\
&=\mathrm{arg}\min\limits_{y\in\mathbb{R}^{k}}\lVert(A+E)\bar{V}_{k}y-(b+\delta_{b})\lVert.\end{split}$

From Theorem 2.1 we know that $\bar{V}_{k}$ is the right orthonormal matrix
generated by the Lanczos biagonalization of $A+E$ with starting vector
$b+\delta_{b}$ in exact arithmetic, which means
$\mathcal{R}(\bar{V}_{k})=\mathcal{K}_{k}((A+E)^{T}(A+E),(A+E)^{T}(b+\delta_{b}))$.
Let $\bar{x}_{k}=\bar{V}_{k}y$. Then $\bar{x}_{k}$ is the $k$-th LSQR solution
to (13). Since $x_{k}=V_{k}y_{k}$, by (10) we have

$\lVert\bar{x}_{k}-x_{k}\lVert\leq\|V_{k}-\bar{V}_{k}\|\|y_{k}\|=O(\|y_{k}\|\mathbf{u})).$

Using $\|y_{k}\|=\|\bar{x}_{k}\|$, we obtain (12).

This result indicates that $x_{k}\approx\bar{x}_{k}$ within $O(\mathbf{u})$.
If we rewrite (13) as

(14) $(A+E)x=(b_{ex}+Ex_{ex})+(e-Ex_{ex}+\delta_{b}),$

then $x_{ex}$ is the exact solution to $(A+E)x_{ex}=b_{ex}+Ex_{ex}$ and
$e-Ex_{ex}+\delta_{b}$ is the noise term. Notice that $\|E\|=O(\mathbf{u})$
and thus the singular values of $A+E$ decrease monotonically without
noticeable gap until they tend to settle at a level of $O(\mathbf{u})$.
Therefore, problem (14) inherits the ill-posedness of (1). Moreover, if
$\|-Ex_{ex}+\delta_{b}\|\ll\|e\|$, then $e-Ex_{ex}+\delta_{b}$ can still be
treated as a Gaussian noise. For a sufficiently small $\mathbf{u}$, we can
hope that the best LSQR regularized solution to (14) has the same accuracy as
that to (1).

In order to compare the accuracy of the best regularized solutions to (1) and
(14), we first review some knowledge concerning the linear ill-posed problem
(1). Since $e\in\mathbb{R}^{m}$ is a Gaussian white noise, all the elements of
$e$ are independent identically distributed (i.i.d.) Gaussian random variables
$\mathcal{N}(0,\sigma^{2})$. Therefore, $\hat{U}^{T}e$ is Gaussian with
$\hat{u}_{i}^{T}e$ are i.i.d. $\mathcal{N}(0,\sigma^{2})$. Thus we have

$m\mathbb{E}[|\hat{u}_{i}^{T}e|^{2}]=\mathbb{E}[\|\hat{U}^{T}e\|^{2}]=\mathbb{E}[\|e\|^{2}]\approx\|e\|^{2}$

and $|\hat{u}_{i}^{T}e|\approx m^{-1/2}\|e\|$; for detailed discussions see
[19, §3.5.1]. The exact right-hand side $b_{ex}$ satisfies the DPC, which can
be written in the following popular simplifying model [17, 19]:

(15) $|\hat{u}_{i}^{T}b_{ex}|=\rho_{0}\sigma_{i}^{1+\beta},\ \ \beta>0,\
i=1,2,\dots,n,$

where $\beta$ is a model parameter that controls the decay rates of
$|\hat{u}_{i}^{T}b_{ex}|$222In [19, §4.6], the corresponding model is
$|u_{i}^{T}b_{ex}|=\sigma_{i}^{1+\beta}$, which does not include the constant
$\rho_{0}$. In fact, Hansen [19, p.68] points out that “while this is, indeed,
a crude model, it reflects the overall behavior often found in real
problems”.. Notice that $\|e\|<\|b_{ex}\|$. The DPC implies that the noisy
coefficients $\hat{u}_{i}^{T}b$ gradually decays in average and are larger
than $|\hat{u}_{i}^{T}e|$ for $i=1,2\dots$ until the noise dominates. Suppose
that the noise in $|\hat{u}_{i}^{T}b|$ starts to dominate at $k=k^{*}+1$,
i.e., $k^{*}$ is the transition point such that [19, §3.5.1]

(16)
$|\hat{u}_{k^{*}}^{T}b|\approx|\hat{u}_{k^{*}}^{T}b_{ex}|>|\hat{u}_{k^{*}}^{T}e|,\
\ |\hat{u}_{k^{*}+1}^{T}b_{ex}|\approx|\hat{u}_{k^{*}+1}^{T}e|.$

Then the DPC of (1) can also be written in the following form:

(17)
$|\hat{u}_{i}^{T}b|=|\hat{u}_{i}^{T}b_{ex}+\hat{u}_{i}^{T}e|\approx\left\\{\begin{array}[]{l}\rho_{0}\sigma_{i}^{1+\beta},\
1\leq i\leq k^{*};\\\ m^{-1/2}\|e\|,\ i\geq k^{*}+1.\end{array}\right.$

The accuracy of the best regularized solution to (1) is closely connected with
the DPC (17), which can be reflected by the effective resolution limit denoted
by $\eta_{res}$. The effective resolution limit of (1) denotes the smallest
coefficient $|\hat{v}_{i}^{T}x_{ex}|$ that can be recovered from the given $A$
and noisy $b$ [17, §4.1], and it is shown in [17, §4.5] that

$\eta_{res}\approx(m^{-1/2}\lVert e\lVert)^{\frac{\beta}{1+\beta}}.$

The accuracy of the best regularized solution $x_{opt}$ is dependent on
$\eta_{res}$, that is, one can only hope that $x_{opt}$ reaches an accuracy
corresponding to the effective resolution limit, which means

(18) $\frac{\|x_{opt}-x_{ex}\|}{\|x_{ex}\|}\geq
C_{1}\varepsilon^{\frac{\beta}{1+\beta}}$

with a moderate constant $C_{1}$, where $\varepsilon=\|e\|/\|b_{ex}\|<1$ is
the noise level; see e.g., [17, 10] for more details.

Now we are ready to analyze the accuracy of computed regularized solutions.
Suppose that the best regularized solution of (14) is $\bar{x}_{k_{0}}$. Let
the $i$-th largest singular values of $A+E$ be $\tilde{\sigma}_{i}$ and the
corresponding left singular vector be $\tilde{u}_{i}$. By the above analysis,
the accuracy of $\bar{x}_{k_{0}}$ will be the same as that of $x_{opt}$ if
333The assertion implicitly uses the order optimal property of the LSQR or its
mathematical equivalent CG algorithm for linear inverse problems [10, §7.3],
which means that the best regularized solution obtained by LSQR can have the
same accuracy as that of $x_{opt}$. This property has not been rigorously
proved for discrete linear ill-posed problems; see [17, §6.3]and [8] for more
discussions.

  * •

the DPC of (14) inherits properties of the DPC of (1), i.e., the DPC of (14)
before the noise dominates should satisfies

(19) $|\tilde{u}_{i}^{T}b|\approx\rho_{0}\tilde{\sigma}_{i}^{1+\beta},\ \
1\leq i\leq k^{*};$

  * •

the effective resolution limit of (14) has the same accuracy as that of (1).

Note that Theorem 3.1 implies that
$\|x_{k_{0}}-\bar{x}_{k_{0}}/\|\bar{x}_{k_{0}}\|=+O(\mathbf{u})$. Therefore,
for the above $\bar{x}_{k_{0}}$, the corresponding $x_{k_{0}}$ will have the
same accuracy as $x_{opt}$ if

(20) $\mathbf{u}\ll C_{1}\varepsilon^{\frac{\beta}{1+\beta}},$

which implies that the best computed regularized solution among $x_{k}$ can
achieve the same accuracy as $x_{opt}$.

In the following we make some analysis about the above assertions.

######  Lemma 3.3.

For the ill-posed problem (13) or its equivalence (14), if the DPC (19) holds
and following two conditions are satisfied: $\mathbf{(1).}$
$\|-Ex_{ex}+\delta_{b}\|\ll\|e\|$; $\mathbf{(2).}$ $\lVert
E\lVert\ll\sigma_{k^{*}}$, then the effective resolution limit of (14) is
approximately equal to that of (1).

######  Proof 3.4.

Since the noise in problem (14) is $e-Ex_{ex}+\delta_{b}$, by Condition
$\mathbf{(1)}$, the noise $e$ dominates, thus we can regard this noise as
Gaussian. Condition $\mathbf{(2)}$ implies that
$\tilde{\sigma}_{i}\approx\sigma_{i}$ until the noise in $|\hat{u}_{i}^{T}b|$
starts to dominates, and thus the errors in $A+E$ starts to dominates at a
point $\bar{k}^{*}>k^{*}$. Therefore, it follows from [17, §4.5], if (19)
holds, that

$\bar{\eta}_{res}\approx(m^{-1/2}\lVert
e-Ex_{ex}+\delta_{b}\lVert)^{\frac{\beta}{1+\beta}}\approx(m^{-1/2}\lVert
e\lVert)^{\frac{\beta}{1+\beta}},$

where $\bar{\eta}_{res}$ is the effective resolution limit of (14).

######  Lemma 3.5.

For the iteration number $k$ not very big, if $\mathbf{u}$ satisfies

(21)
$\mathbf{u}\ll\min\\{\varepsilon,(m^{-1/2}\varepsilon)^{\frac{1}{1+\beta}}\\},$

then the relation (20) and Conditions $\mathbf{(1)}$, $\mathbf{(2)}$ hold.

######  Proof 3.6.

Condition $\mathbf{(1)}$ holds if $\|Ex_{ex}\|\ll\|e\|$ and
$\|\delta_{b}\|\ll\|e\|$. By (11), these two equalities can be satisfied if
$\mathbf{u}\ll\|e\|/\|b\|$ and $\|E\|\ll\|e\|/\|x_{ex}\|$. Since

$\dfrac{\|e\|}{\|b\|}\geq\dfrac{\|e\|}{\|b_{ex}\|+\|e\|}=\dfrac{\varepsilon}{1+\varepsilon}>\dfrac{\varepsilon}{2},$

we have $\mathbf{u}\ll\|e\|/\|b\|$ if $\mathbf{u}\ll\varepsilon$. Note that

$\dfrac{\|e\|}{\|x_{ex}\|}\leq\dfrac{\|A\|\|e\|}{\|b_{ex}\|}=\varepsilon\|A\|,$

and the value of $\|e\|/\|x_{ex}\|$ should not deviate too far from
$\varepsilon\|A\|$ since $\|x_{ex}\|$ is usually a moderate quantity. Using
$\lVert E\lVert=O(c(n,k)\|A\|\mathbf{u})$, we have $\|E\|\ll\|e\|/\|x_{ex}\|$
if $\mathbf{u}\ll\varepsilon$ since $c(n,k)$ is a moderate quantity.

By (15) and (16), we have
$\rho_{0}\sigma_{k^{*}+1}^{1+\beta}\approx|\hat{u}_{k^{*}+1}^{T}e|\approx
m^{-1/2}\|e\|$, which implies that Condition $\mathbf{(2)}$ will hold if

(22)
$\|E\|\ll\sigma_{k^{*}+1}\approx(m^{-1/2}\rho_{0}^{-1}\|e\|)^{\frac{1}{1+\beta}}.$

By (15), we have $\rho_{0}=|\hat{u}_{1}^{T}b_{ex}|/\sigma_{1}^{1+\beta}$, and
thus

(23)
$(\rho_{0}^{-1}\|e\|)^{\frac{1}{1+\beta}}=\sigma_{1}\Big{(}\dfrac{\|e\|}{|\hat{u}_{1}^{T}b_{ex}|}\Big{)}^{\frac{1}{1+\beta}}\geq\|A\|\Big{(}\dfrac{\|e\|}{\|b_{ex}\|}\Big{)}^{\frac{1}{1+\beta}}=\varepsilon^{\frac{1}{1+\beta}}\|A\|.$

Therefore Condition $\mathbf{(2)}$ will hold if
$\|E\|\ll(m^{-1/2}\varepsilon)^{\frac{1}{1+\beta}}\|A\|$, which can be
satisfied if $\mathbf{u}\ll(m^{-1/2}\varepsilon)^{\frac{1}{1+\beta}}$. By the
above derivations, one can check that (20) and Condition $\mathbf{(1)}$,
$\mathbf{(2)}$ hold if $\mathbf{u}$ satisfies (21).

In order to analyze the DPC (19), we adopt the following popular model to
describe the decay rates of $\sigma_{i}$ for different types of ill-posedness
[17]:

(24) $\sigma_{i}=\left\\{\begin{array}[]{l}\zeta\rho^{-i},\ \rho>1\ \ \
\mathrm{severely\ ill\mbox{-}posed};\\\ \zeta i^{-\alpha},\ \alpha>1\ \ \
\mathrm{moderately\ ill\mbox{-}posed};\\\ \zeta i^{-\alpha},\ 1/2<\alpha\leq
1\ \ \ \mathrm{mildly\ ill\mbox{-}posed}.\end{array}\right.$

######  Remark 3.7.

Notice that the model (24) means that all the singular values of $A$ are
simple. For the case that $A$ has multiple singular values, the model can be
rewritten by the following modification. First rewrite the SVD of $A$ as

$A=\widehat{U}\begin{pmatrix}\Sigma\\\
\mathbf{0}\end{pmatrix}\widehat{V}^{T},$

where
$\widehat{U}=(\widehat{U}_{1},\dots,\widehat{U}_{r},\widehat{U}_{\perp})$ with
$\widehat{U}_{i}\in\mathbb{R}^{m\times l_{i}}$ and
$\widehat{V}=(\widehat{V}_{1},\dots,\widehat{V}_{r})$ with
$\widehat{V}_{i}\in\mathbb{R}^{n\times l_{i}}$ are column orthonormal,
$\Sigma=\mbox{diag}(\hat{\sigma}_{1}I_{l_{1}},\dots,\hat{\sigma}_{r}I_{l_{r}})$
with the $r$ distinct singular values
$\hat{\sigma}_{1}>\hat{\sigma}_{2}>\dots>\hat{\sigma}_{r}>0$, each
$\hat{\sigma}_{i}$ is $l_{i}$ multiple and $l_{1}+l_{2}+\cdots+l_{r}=n$. Then
the decay rate of $\hat{\sigma}_{i}$ can be written in the same form as (24).
In this case, the DPC of (1) becomes

(25) $\|\widehat{U}_{i}^{T}b_{ex}\|=\rho_{0}\hat{\sigma}_{i}^{1+\beta},\ \
i=1,2,\dots,r,$

which states that, on average the (generalized) Fourier coefficients
$\|\widehat{U}_{i}^{T}b_{ex}\|$ decay faster than $\hat{\sigma}_{i}$.

Using this model, we can give a sufficient condition under which the relation
(19) holds. For notational simplicity, we also write $\hat{\sigma}_{i}$ as
$\sigma_{i}$ without causing confusions.

######  Lemma 3.8.

Suppose that the iteration number $k$ is not very big and
$\sigma_{i}-\sigma_{i+1}\gg\lVert E\lVert$ for $1\leq i\leq k^{*}$. If (21)
holds and

(26)
$\mathbf{u}\ll(m^{-1/2}\varepsilon)^{\frac{2+\beta}{1+\beta}}\big{(}\frac{\sigma_{k^{*}}}{\sigma_{k^{*}+1}}-1\big{)},$

then the DPC of (14) satisfies (19).

######  Proof 3.9.

There two cases needed to be proved.

Case 1. $A$ has single singular values. Write the $i$-th left singular vector
of $A+E$ as $\tilde{u}_{i}=u_{i}+\delta_{u_{i}}$ where $\delta_{u_{i}}$ is an
error vector. Since (21) holds, we have $\tilde{\sigma}_{i}\approx\sigma_{i}$
for $1\leq i\leq k^{*}$. Note (19) implies
$|(u_{i}+\delta_{u_{i}})^{T}b|\approx\rho_{0}\tilde{\sigma}_{i}^{1+\beta}\approx\rho_{0}\sigma_{i}^{1+\beta}$
for $1\leq i\leq k^{*}$, which can be satisfied if

(27) $|\delta_{u_{i}}^{T}b|\ll\rho_{0}\sigma_{i}^{1+\beta},\ \ 1\leq i\leq
k^{*}.$

By the perturbation theorem of singular vectors [4, Theorem 1.2.8], we have
the perturbation bound

$|\sin\theta(u_{i},\tilde{u}_{i})|\leq\frac{\lVert
E\lVert}{\sigma_{i}-\sigma_{i+1}-\lVert E\lVert}\approx\frac{\lVert
E\lVert}{\sigma_{i}-c_{i+1}},\ \ 1\leq i\leq k^{*}$

under the assumption that $\sigma_{i}-\sigma_{i+1}\gg\lVert E\lVert$ for
$1\leq i\leq k^{*}$, where $\theta(u_{i},\tilde{u}_{i})$ is the angle between
$u_{i}$ and $\tilde{u}_{i}$. Thus we have

$\lVert\delta_{u_{i}}\lVert=2|\sin(\theta(u_{i},\tilde{u}_{i})/2)|\leq|\sin\theta(u_{i},\tilde{u}_{i})|\lesssim\frac{\lVert
E\lVert}{\sigma_{i}-\sigma_{i+1}}.$

Therefore, (27) can be satisfied if

$\frac{\lVert E\lVert\lVert
b\lVert}{\sigma_{i}-\sigma_{i+1}}\ll\rho_{0}\sigma_{i}^{1+\beta},\ \ 1\leq
i\leq k^{*},$

which is equivalent to

(28) $\lVert
E\lVert\ll\frac{\rho_{0}\sigma_{i}^{2+\beta}\big{(}1-\frac{\sigma_{i+1}}{\sigma_{i}}\big{)}}{\lVert
b\lVert},\ \ 1\leq i\leq k^{*}.$

Using the expression of $\sigma_{k^{*}+1}$ in (22) and the model (24), the
minimum of the right-hand term of the above inequality is achieved at
$i=k^{*}$, which is

$\displaystyle\frac{\rho_{0}\sigma_{k^{*}}^{2+\beta}\big{(}1-\frac{\sigma_{k^{*}+1}}{\sigma_{k^{*}}}\big{)}}{\lVert
b\lVert}$
$\displaystyle=\frac{\rho_{0}\sigma_{k^{*}+1}^{2+\beta}\big{(}\frac{\sigma_{k^{*}}}{\sigma_{k^{*}+1}}\big{)}^{2+\beta}\big{(}1-\frac{\sigma_{k^{*}+1}}{\sigma_{k^{*}}}\big{)}}{\lVert
b\lVert}$
$\displaystyle\approx\frac{\rho_{0}(m^{-1/2}\rho_{0}^{-1}\|e\|)^{\frac{2+\beta}{1+\beta}}\big{(}\frac{\sigma_{k^{*}}}{\sigma_{k^{*}+1}}\big{)}^{1+\beta}\big{(}\frac{\sigma_{k^{*}}}{\sigma_{k^{*}+1}}-1\big{)}}{\lVert
b_{ex}\lVert}$
$\displaystyle\geq\frac{(m^{-1/2}\|e\|)^{\frac{2+\beta}{1+\beta}}\big{(}\frac{\sigma_{k^{*}}}{\sigma_{k^{*}+1}}-1\big{)}}{\lVert
b_{ex}\lVert\rho_{0}^{\frac{1}{1+\beta}}}.$

By (23), we have

$\displaystyle\frac{(m^{-1/2}\|e\|)^{\frac{2+\beta}{1+\beta}}}{\lVert
b_{ex}\lVert\rho_{0}^{\frac{1}{1+\beta}}}$
$\displaystyle=(m^{-1/2})^{\frac{2+\beta}{1+\beta}}(\rho_{0}^{-1}\|e\|)^{\frac{1}{1+\beta}}\frac{\|e\|}{\lVert
b_{ex}\lVert}$
$\displaystyle\geq(m^{-1/2})^{\frac{2+\beta}{1+\beta}}\varepsilon^{\frac{1}{1+\beta}}\|A\|\varepsilon=(m^{-1/2}\varepsilon)^{\frac{2+\beta}{1+\beta}}\|A\|.$

Therefore, by (27) and (28) and using $\lVert
E\lVert=O(c(n,k)\|A\|\mathbf{u})$, we finally obtain the result.

Case 2. $A$ has multiple singular values. Write the SVD of $A+E$ in a similar
form as that of $A$, such that the left singular vectors can be written as
$\tilde{U}=(\tilde{U}_{1},\dots,\tilde{U}_{r},\tilde{U}_{\perp})$. By the
perturbation theorem of invariant singular subspaces [38], we have

$\|\sin\Theta(\widehat{U}_{i},\tilde{U}_{i})\|\leq\frac{\lVert
E\lVert}{\hat{\sigma}_{i}-\hat{\sigma}_{i+1}-\lVert
E\lVert}\approx\frac{\lVert E\lVert}{\hat{\sigma}_{i}-\hat{\sigma}_{i+1}}$

where
$\|\sin\Theta(\widehat{U}_{i},\tilde{U}_{i})\|=\lVert\widehat{U}_{i}\widehat{U}_{i}^{T}-\tilde{U}_{i}\tilde{U}_{i}^{T}\lVert$
is the angle measure between subspaces spanned by $\widehat{U}_{i}$ and
$\tilde{U}_{i}$ [13, §2.5]. Notice that

$\displaystyle|\lVert\tilde{U}_{i}^{T}b\lVert-\lVert\widehat{U}_{i}^{T}b\lVert|$
$\displaystyle=|\lVert\tilde{U}_{i}\tilde{U}_{i}^{T}b\lVert-\lVert\widehat{U}_{i}\widehat{U}_{i}^{T}b\lVert|\leq\lVert(\tilde{U}_{i}\tilde{U}_{i}^{T}-\widehat{U}_{i}\widehat{U}_{i}^{T})b\lVert$
$\displaystyle\leq\lVert
b\lVert\|\sin\Theta(\widehat{U}_{i},\tilde{U}_{i})\|,$

where $|\lVert\tilde{U}_{i}^{T}b\lVert-\lVert\widehat{U}_{i}^{T}b\lVert|$ is
the corresponding version of the left-hand term of (27). Using the same
approach as that for analyzing (27), we can obtain the result for the multiple
singular values case.

Using model (24), the minimum of $\sigma_{i}-\sigma_{i+1}$ for $1\leq i\leq k$
is achieved at $i=k^{*}$. Thus for $1\leq i\leq k^{*}$, we have

$\sigma_{i}-\sigma_{i+1}=\sigma_{i+1}\big{(}\frac{\sigma_{i}}{\sigma_{i+1}}-1\big{)}\geq\left\\{\begin{array}[]{l}\sigma_{k^{*}+1}(\rho-1)\
\ \ \mathrm{severely\ ill\mbox{-}posed};\\\
\sigma_{k^{*}+1}[(\frac{k^{*}+1}{k^{*}})^{\alpha}-1]\ \ \
\mathrm{moderately/mildly\ ill\mbox{-}posed}.\end{array}\right.$

By (22) and (23), we have
$\sigma_{k^{*}+1}\gtrsim(m^{-1/2}\varepsilon)^{\frac{1}{1+\beta}}\|A\|$. Using
these two inequalities, one can check that if (21) and (26) hold, then the
assumption that $\sigma_{i}-\sigma_{i+1}\gg\lVert E\lVert$ for $1\leq i\leq
k^{*}$ can be satisfied.

Note that the semi-convergence point $k_{0}$ of LSQR for (14) is usually not
big and thus $\lVert E\lVert=O(c(n,k_{0})\|A\|\mathbf{u})$ with $c(n,k_{0})$ a
moderate quantity. By Lemma 3.5 and Lemma 3.8, if (21) and (26) hold, then
relation (20) holds and the effective resolution limit of (14) is
approximately equal to that of (1). Therefore $\bar{x}_{k_{0}}$ as well as
$x_{k_{0}}$ will have the same accuracy as $x_{opt}$, which implies that the
best computed regularized solution among $x_{k}$ can have the same accuracy as
$x_{opt}$. We finally obtain the following result.

######  Theorem 3.10.

For the LSQR with Lanczos bidiagonalization (with full reorthogonalization)
implemented in finite precision arithmetic, if the roundoff unit $\mathbf{u}$
satisfies

(29) $\mathbf{u}\ll\varrho(m^{-1/2}\varepsilon)^{\frac{2+\beta}{1+\beta}}$

where

(30) $\varrho=\left\\{\begin{array}[]{l}\min\\{1,\rho-1\\}\ \ \
\mathrm{severely\ ill\mbox{-}posed};\\\
\min\\{1,(\frac{k^{*}+1}{k^{*}})^{\alpha}-1\\}\ \ \ \mathrm{moderately/mildly\
ill\mbox{-}posed},\end{array}\right.$

then the best computed regularized solution to (1) can achieve the same
accuracy as the best regularized LSQR solution to (1) obtained in exact
arithmetic.

In Theorem 3.10, the parameters $\alpha$, $\beta$ and $\rho$ are unknown in
practical computations. In fact, these parameters are ideal for simplifying
singular value decaying and DPC models, and they are closely related to
properties of a given ill-posed problem. However, it is instructive from them
to get insight into a practical choice of $\mathbf{u}$. For severely ill-posed
problems, $\rho-1$ is usually a constant not very small, while for
moderately/mildly ill-posed problems, if $\varepsilon$ is very small and
$\alpha$ is not big that means the singular values of $A$ decaying very
slowly, then $k^{*}$ will be big and thus
$(\frac{k^{*}+1}{k^{*}})^{\alpha}-1\approx\alpha/k^{*}$ will be very small.
Therefore, for moderately/mildly ill-posed problems, if $\varepsilon$ is very
small and $\alpha$ is not big, then (29) may give a too small upper bound on
$\mathbf{u}$.

Theorem 3.10 implies that for noisy level $\varepsilon$ not very small, we can
exploit lower precision for constructing Lanczos vectors in the LSQR for
solving (1) without loss of any accuracy of final regularized solutions. We
will use numerical examples to show that single precision is enough for the
three types of linear ill-posed problem. It is important to stress a special
practical case that $k^{*}=n$ for a too small $\varepsilon$ and $\alpha$,
which may be encountered in some image deblurring problems. In this case the
noise amplification is tolerable even without regularization, which makes
$x_{nai}$ a good approximation to $x_{ex}$, and the LSQR solves (1) in their
standard manners as if they solved an ordinary other than ill-posed problem.
Thus our result can not be applied to this case.

##  4 Updating $x_{k}$ using lower precision

In this section, we discuss how to use lower precision for updating $x_{k}$
step by step. Suppose that the $k$-step Lanczos bidiagonalization is
implemented using computing precision chosen as in Theorem 3.10. We first
review the procedure for updating $x_{k}$ from $x_{0}=\mathbf{0}$ proposed in
[32]. First, the QR factorization

$Q_{k}\begin{pmatrix}B_{k}&\beta_{1}e_{1}^{(k+1)}\end{pmatrix}=\begin{pmatrix}R_{k}&f_{k}\\\
&\bar{\phi}_{k+1}\end{pmatrix}=\left(\begin{array}[]{ccccc;{2pt/2pt}c}\rho_{1}&\theta_{2}&&&&\phi_1\\\
&\rho_{2}&\theta_{3}&&&\phi_2\\\ &&\ddots&\ddots&&\vdots\\\
&&&\rho_{k-1}&\theta_{k}&\phi_{k-1}\\\ &&&&\rho_{k}&\phi_k\\\
\hdashline[2pt/2pt]&&&&&\bar{\phi}_{k+1}\end{array}\right)$

is performed using a series of Givens rotations, where at the $i$-th step the
Givens rotation is chosen to zero out $\beta_{i+1}$:

$\begin{pmatrix}c_{i}&s_{i}\\\
s_{i}&-c_{i}\end{pmatrix}\begin{pmatrix}\bar{\rho}_{i}&0&\bar{\phi}_{i}\\\
\beta_{i+1}&\alpha_{i+1}&0\end{pmatrix}\begin{pmatrix}\rho_{i}&\theta_{i+1}&\phi_{i}\\\
0&\bar{\rho}_{i+1}&\bar{\phi}_{i+1}\end{pmatrix},$

and the orthogonal matrix $Q_{k}$ is the product of these Givens rotation
matrices. Since

(31)
$\|B_{k}y-\beta_{1}e_{1}^{(k+1)}\|^{2}=\Big{\|}Q_{k}\begin{pmatrix}B_{k}&\beta_{1}e_{1}^{(k+1)}\end{pmatrix}\begin{pmatrix}y\\\
-1\end{pmatrix}\Big{\|}^{2}=\|R_{k}y-f_{k}\|^{2}+|\bar{\phi}_{k+1}|^{2},$

the solution to $\min_{y\in\mathbb{R}^{k}}\|B_{k}y-\beta_{1}e_{1}^{(k)}\|$ is
$y_{k}=R_{k}^{-1}f_{k}$. Factorize $R_{k}$ as

$R_{k}=D_{k}\widetilde{R}_{k},\ \ D_{k}=\begin{pmatrix}\rho_{1}&&&\\\
&\rho_{2}&&\\\ &&\ddots&\\\ &&&\rho_{k}\end{pmatrix},\
\widetilde{R}_{k}=\begin{pmatrix}1&\theta_{2}/\rho_{1}&&\\\
&1&\theta_{3}/\rho_{2}&\\\ &&\ddots&\theta_{k}/\rho_{k-1}\\\
&&&1\end{pmatrix},$

then we get

(32)
$x_{k}=V_{k}y_{k}=V_{k}R_{k}^{-1}f_{k}=(V_{k}\widetilde{R}_{k}^{-1})(D_{k}^{-1}f_{k}).$

Let $W_{k}=V_{k}\widetilde{R}_{k}^{-1}=(w_{1},\dots,w_{k})$. By using back
substitution for solving $W_{k}\widetilde{R}_{k}=V_{k}$ we obtain the updating
procedure for $x_{i}$ and $w_{i}$:

(33) $x_{i}=x_{i-1}+(\phi_{i}/\rho_{i})w_{i},\ \ \
w_{i+1}=v_{i+1}-(\theta_{i+1}/\rho_{i})w_{i},$

which is described in Algorithm 2.

Algorithm 2 Updating procedure

1: Let $x_{0}=\mathbf{0},\ w_{1}=v_{1},\ \bar{\phi}_{1}=\beta_{1},\
\bar{\rho}_{1}=\alpha_{1}$

2: for $i=1,2,\ldots,k,$ do

3: $\rho_{i}=(\bar{\rho}_{i}^{2}+\beta_{i+1}^{2})^{1/2}$

4: $c_{i}=\bar{\rho}_{i}/\rho_{i},\ s_{i}=\beta_{i+1}/\rho_{i}$

5: $\theta_{i+1}=s_{i}\alpha_{i+1},\ \bar{\rho}_{i+1}=-c_{i}\alpha_{i+1}$

6: $\phi_{i}=c_{i}\bar{\phi}_{i},\ \bar{\phi}_{i+1}=s_{i}\bar{\phi}_{i}$

7: $x_{i}=x_{i-1}+(\phi_{i}/\rho_{i})w_{i}$

8: $w_{i+1}=v_{i+1}-(\theta_{i+1}/\rho_{i})w_{i}$

9: end for

In Algorithm 2 the Givens QR factorization (Line 3 – Line 6) can be performed
very quickly since only floating point arithmetic of scalars is involved. In
contrast, the updating of $x_{i}$ and $w_{i+1}$ involves vector computations.
Therefore, in practical computation, we can use double precision to perform
the $QR$ factorization, while the updating of $x_{i}$ and $w_{i+1}$ are
performed using lower precision. In order to simplify the theoretical analysis
for this implementation of the updating procedure, we assume that the QR
factorization is performed in exact arithmetic, while $x_{i}$ and $w_{i+1}$
are computed in finite precision arithmetic with roundoff unit
$\mathbf{\bar{u}}$.

We use the following model [13, §2.7.3] for the floating point arithmetic:

(34) $\displaystyle\mbox{fl}(a\ \mbox{op}\ b)=(a\ \mbox{op}\ b)(1+\epsilon),\
\ \ |\epsilon|\leq\mathbf{\bar{u}},\ \ \mbox{op}=+,-,*,/.$

Under this model, we have the following rounding error results for matrix and
vector computations [13, §2.7.8]:

(35) $\displaystyle\mbox{fl}(u+\alpha v)=u+\alpha v+w,\ \ \
|w|\leq(|u|+2|\alpha v|)\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}),$ (36)
$\displaystyle\mbox{fl}(AB)=AB+X,\ \ \ |X|\leq
n|A||B|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}).$

where $u,v$ are vectors, $\alpha$ is a scalar, and $A,B$ are two matrices of
orders $m\times n$ and $n\times l$, respectively. In (35) and (36), the
notation $|\cdot|$ is used to denote the absolute value of a matrix or vector
and $``\leq"$ means the relation $``\leq"$ holds componentwise. We remark that
(36) applies to both dot-product and outer-product based procedures for matrix
multiplications; see [13, §1.1, §2.7] for these two types of computation of
matrix multiplications and the corresponding rounding error results.

Denote the computed version of $w_{i}$ by $\hat{w}_{i}$. Using the above
model, the computation of $\hat{w}_{i}$ in finite precision arithmetic can be
formed as:

$\left\\{\begin{array}[]{l}\hat{w}_{1}=v_{1}+\delta_{1}^{w}\\\
\hat{w}_{2}=v_{2}-(\theta_{2}/\rho_{1})\hat{w}_{1}+\delta_{2}^{w}\\\ \ \ \ \ \
\vdots\\\
\hat{w}_{k}=v_{k}-(\theta_{k}/\rho_{k-1})\hat{w}_{k-1}+\delta_{k}^{w}\end{array}\right.$

where

$\|\delta_{1}^{w}\|\leq\|v_{1}\|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})=\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})$

and

$\|\delta_{i}^{w}\|\leq(\|v_{i}\|+2\|(\theta_{i}/\rho_{i-1})\hat{w}_{i-1}\|)\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})=(1+2\|(\theta_{i}/\rho_{i-1})\hat{w}_{i-1}\|)\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}).$

Let $h_{i}=\hat{w}_{i}-w_{i}$. Then we have

$h_{i+1}=v_{i+1}-(\theta_{i+1}/\rho_{i})\hat{w}_{i}+\delta_{i}^{w}-(v_{i+1}-(\theta_{i+1}/\rho_{i})w_{i})=-(\theta_{i+1}/\rho_{i})h_{i}+\delta_{i}^{w},$

which leads to

(37) $H_{k}\widetilde{R}_{k}=\Delta_{k}^{w}$

with $H_{k}=(h_{1},\dots,h_{k})$ and
$\Delta_{k}^{w}=(\delta_{1}^{w},\dots,\Delta_{k}^{w})$. Therefore, we have

$\displaystyle\|\Delta_{k}^{w}\|$ $\displaystyle\leq\sqrt{k}\max_{1\leq i\leq
k}\|\delta_{i}^{w}\|\leq\sqrt{k}[1+2\max_{1\leq i\leq
k}(\theta_{i+1}/\rho_{i})\|w_{i}+h_{i}\|]\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})$
$\displaystyle\leq\sqrt{k}[1+2\|\widetilde{R}_{k}\|(\|W_{k}\|+\|\Delta_{k}^{w}\|\|\widetilde{R}_{k}^{-1}\|)]\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}),$

where we have used $|\theta_{i+1}/\rho_{i}|\leq\|\widetilde{R}_{k}\|$ and
$\|h_{i}\|\leq\|H_{k}\|\leq\|\Delta_{k}^{w}\|\|\widetilde{R}_{k}^{-1}\|$. This
inequality leads to

$(1-2\sqrt{k}\kappa(\widetilde{R}_{k})\mathbf{\bar{u}})\|\Delta_{k}^{w}\|\leq\sqrt{k}(1+2\|\widetilde{R}_{k}\|\|W_{k}\|)\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}),$

where
$\kappa(\widetilde{R}_{k})=\|\widetilde{R}_{k}\|\|\widetilde{R}_{k}^{-1}\|$ is
the condition number of $\widetilde{R}_{k}$. Notice that
$\|W_{k}\|\leq\|V_{k}\|\|\widetilde{R}_{k}^{-1}\|=\|\widetilde{R}_{k}^{-1}\|+O(\mathbf{u})$
where we have used $\|V_{k}-\bar{V}_{k}\|=O(\mathbf{u})$ and
$\|\bar{V}_{k}\|=1$ by (10). We obtain the upper bound on
$\|\Delta_{k}^{w}\|$:

(38)
$\|\Delta_{k}^{w}\|\leq\sqrt{k}[1+2(1+\sqrt{k})\kappa(\widetilde{R}_{k})]\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}+\mathbf{\bar{u}}\mathbf{u}).$

Now we can analyze the accuracy of the computed $x_{k}$ updated using
computing precision with roundoff unit $\mathbf{\bar{u}}$.

######  Theorem 4.1.

Suppose the Lanczos bidiagonalization (with full reorthogonalization) in LSQR
is implemented in finite precision arithmetic with roundoff unit $\mathbf{u}$
and $x_{k}$ is the exact solution to (6). Denote the practical updated $x_{k}$
by $\hat{x}_{k}$. We have

(39)
$\dfrac{\|\hat{x}_{k}-x_{k}\|}{\|x_{k}\|}\leq\sqrt{k}[1+(2+2\sqrt{k}+k)\kappa(\widetilde{R}_{k})]\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}+\mathbf{\bar{u}}\mathbf{u}).$

######  Proof 4.2.

Notice from (32) and (33) that the formation of $\hat{x}_{k}$ is the matrix
multiplication between $\widehat{W}_{k}=(\hat{w}_{1},\dots,\hat{w}_{k})$ and
$D_{k}^{-1}f_{k}$ by the outer-product based procedure. Using model (36) we
have

(40)
$\hat{x}_{k}-x_{k}=\widehat{W}_{k}D_{k}^{-1}f_{k}+\Delta_{k}^{x}-W_{k}D_{k}^{-1}f_{k},$

where $\hat{x}_{k}=\widehat{W}_{k}(D_{k}^{-1}f_{k})+\Delta_{k}^{x}$ with
$|\Delta_{k}^{x}|\leq
k|\widehat{W}_{k}||D_{k}^{-1}f_{k}|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})$.
Therefore, we have

$\displaystyle\|\Delta_{k}^{x}\|$ $\displaystyle\leq\||\Delta_{k}^{x}|\|\leq
k\||\widehat{W}_{k}|\|_{F}\||D_{k}^{-1}f_{k}|\|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})$
$\displaystyle=k\|\widehat{W}_{k}\|_{F}\|D_{k}^{-1}f_{k}\|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})$
$\displaystyle\leq
k^{3/2}\|\widehat{W}_{k}\|\|D_{k}^{-1}f_{k}\|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}),$

where $\|\cdot\|_{F}$ is the Frobenius norm of a matrix. By (37), we have
$\widehat{W}_{k}-W_{k}=H_{k}=\Delta_{k}^{w}\widetilde{R}_{k}^{-1}$, and thus

$\|\widehat{W}_{k}\|\leq\|W_{k}\|+\|H_{k}\|\leq\|W_{k}\|+\|\Delta_{k}^{w}\|\|\widetilde{R}_{k}^{-1}\|.$

Substituting it into the inequality about $\|\Delta_{k}^{x}\|$ and noticing
$\widetilde{R}_{k}^{-1}D_{k}^{-1}f_{k}=R_{k}^{-1}f_{k}=y_{k}$, we obtain

$\displaystyle\|\Delta_{k}^{x}\|$ $\displaystyle\leq
k^{3/2}(\|W_{k}\|+\|\Delta_{k}^{w}\|\|\widetilde{R}_{k}^{-1}\|)\|\widetilde{R}_{k}(\widetilde{R}_{k}^{-1}D_{k}^{-1}f_{k})\|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2})$
$\displaystyle\leq
k^{3/2}\kappa(\widetilde{R}_{k})\|y_{k}\|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}+\mathbf{\bar{u}}\mathbf{u}),$

where we have used $\|W_{k}\|\leq\|\widetilde{R}_{k}^{-1}\|+O(\mathbf{u})$.
Using $\widehat{W}_{k}-W_{k}=\Delta_{k}^{w}\widetilde{R}_{k}^{-1}$ again, we
get

$\displaystyle\|\widehat{W}_{k}D_{k}^{-1}f_{k}-W_{k}D_{k}^{-1}f_{k}\|=\|(\widehat{W}_{k}-W_{k})D_{k}^{-1}f_{k}\|=\|\Delta_{k}^{w}(\widetilde{R}_{k}^{-1}D_{k}^{-1}f_{k})\|\leq\|\Delta_{k}^{w}\|\|y_{k}\|.$

By (40) and combining with (38), we obtain

$\displaystyle\|\hat{x}_{k}-x_{k}\|$
$\displaystyle\leq\|\Delta_{k}^{x}\|+\|\widehat{W}_{k}D_{k}^{-1}f_{k}-W_{k}D_{k}^{-1}f_{k}\|$
$\displaystyle\leq\sqrt{k}[1+(2+2\sqrt{k}+k)\kappa(\widetilde{R}_{k})]\|y_{k}\|\mathbf{\bar{u}}+O(\mathbf{\bar{u}}^{2}+\mathbf{\bar{u}}\mathbf{u}).$

Noticing that $x_{k}$ is the exact solution to (6), where the $B_{k}$ and
Lanczos vectors are computed in finite precision arithmetic with roundoff unit
$\mathbf{u}$. By using the relation (12) and
$\|y_{k}\|=\|\bar{V}_{k}y_{k}\|=\|\bar{x}_{k}\|$, we finally obtain (39).

Note that $\widetilde{R}_{k}$ is obtained by scaling $R_{k}$ using the
diagonal of it, and the diagonal scaling step can often dramatically reduces
the condition number of $R_{k}$. In fact, we will show in the numerical
experiments section that $\kappa(\widetilde{R}_{k})$ is a moderate value even
for an iteration $k$ bigger than the semi-convergence point. By Theorem 4.1,
when $\mathbf{u}$ has been chosen such that $x_{k}$ can achieve the same
accuracy of the best LSQR regularized solution to (1) obtained in exact
arithmetic, in order to make the practical updated $\hat{x}_{k}$ can also
achieve the same accuracy, $\mathbf{\bar{u}}$ should be chosen such that the
upper bound in (39) is much smaller than $\|x_{opt}-x_{ex}\|/\|x_{ex}\|$. For
a not very small noise level, the single precision roundoff unit
$\mathbf{\bar{u}}$ is enough. This ensures that we can use lower precision for
updating $x_{k}$, which is more efficient than using double precision.

Finally we discuss methods for estimating the optimal early stopping
iteration. By Theorem 3.1 and (31) we have

$\bar{\phi}_{k+1}=\|B_{k}y_{k}-\beta_{1}e_{1}^{(k+1)}\|=\lVert(A+E)x_{k}-(b+\delta_{b})\lVert.$

For the proper choice of $\mathbf{u}$, the noise norm of (13) is
$\|e-Ex_{ex}+\delta_{b}\|\approx\|e\|$ since
$\|-Ex_{ex}+\delta_{b}\|\ll\|e\|$. Notice that the computed $x_{k}$ in lower
precision is actually the solution to (13) and $\bar{\phi}_{k+1}$ is computed
before practically updating $\hat{x}_{k}$. Therefore, the discrepancy
principle corresponding to (13) can be written as

$\lVert(A+E)x_{k}-(b+\delta_{b})\lVert\lesssim\tau\|e\|$

with $\tau>1$ slightly, and we should stop iteration at the first $k$
satisfying

(41) $\bar{\phi}_{k+1}=\|B_{k}y_{k}-\beta_{1}e_{1}^{(k+1)}\|\leq\tau\|e\|,$

and use this $k$ as the estimate of $k_{0}$, where $\bar{\phi}_{k+1}$ can be
efficiently computed by using Algorithm 2. Numerical experiments will show
that this estimate is almost the same as that obtained by the discrepancy
principle for LSQR in double precision arithmetic. The discrepancy principle
method usually suffers from under-estimate $k_{0}$ and thus the solution is
over-regularized.

Another approach is the L-curve criterion, which does not need $\|e\|$ in
advance. The motivation is that one can plot
$\left(\log\|Ax_{k}-b\|,\log\|x_{k}\|\right)$ in the shape of an L-curve, and
the corner of the curve is a good estimate of the optimal early stopping
iteration. The L-curve for (13) is

$\left(\log\|(A+E)x_{k}-(b+\delta_{b})\|,\log\|x_{k}\|\right),$

which is just

(42) $\left(\log\bar{\phi}_{k+1},\log\|x_{k}\|\right),$

where the norm of $x_{k}$ should be computed at each iteration. If we use
lower precision for updating $x_{k}$, we should compute the norm of
$\hat{x}_{k}$ instead of $x_{k}$, and this may make little difference with the
estimate of $k_{0}$. Numerical experiments will show that this estimate is
almost the same as that obtained by the L-curve criterion for LSQR in double
precision arithmetic.

##  5 Numerical Experiments

In this section, we present some numerical experiments to illustrate
performance of LSQR for discrete linear ill-posed problems implemented using
different computing precisions. We use “d” to denote the algorithm implemented
using double precision, and use “s+d” and “s+s” to denote the algorithms that
use single precision for the Lanczos bidiagonalization while use double and
single precisions for updating $x_{k}$, respectively. For these different
implementations, we compare accuracy of the regularized solutions by using the
relative reconstruction error

$RE(k)=\frac{\|x_{k}-x_{ex}\|}{\|x_{ex}\|}$

to plot semi-convergence curves, where $x_{k}$ denote the computed solutions
produced by the three implementations.

We choose some one dimensional (1-D) problems from the regularization toolbox
[18], and two dimensional (2-D) image deblurring problems from [11]. The
description of all test examples is listed in Table 5.1. All experiments are
performed in MATLAB R2019b, where the roundoff units for double and single
precision arithmetic are $2^{-53}\approx 1.11\times 10^{-16}$ and
$2^{-24}\approx 5.96\times 10^{-8}$, respectively. The codes are available at
https://github.com/Machealb/Lower_precision_solver.

Table 5.1: The description of test problems Problem | $m\times n$ | Ill-posedness | Description
---|---|---|---
shaw | $1000\times 1000$ | severe | 1-D image restoration model
deriv2 | $1000\times 1000$ | moderate | Computation of second derivative
gravity | $2000\times 2000$ | severe | 1-D gravity surveying problem
heat | $2000\times 2000$ | moderate | Inverse heat equation
PRblurspeckle | $16384\times 16384$ | mild | 2-D image deblurring problem
PRblurdefocus | $65536\times 65536$ | mild | 2-D image deblurring problem

###  5.1 One dimensional case

For one dimensional problems shaw, deriv2, gravity and heat, we use the codes
from [18] to generate $A$, $x_{ex}$ and $b_{ex}=Ax_{ex}$, and then add a white
Gaussian noise $e$ with a prescribed noise level
$\varepsilon=\|e\|/\|b_{ex}\|$ to $b_{ex}$ and form the noisy $b=b_{ex}+e$.

(a) shaw

(b) deriv2

(c) gravity

(d) heat

Figure 5.1: Semi-convergence curves for LSQR implemented using different
computing precisions, $\varepsilon=10^{-3}$.

(a) shaw

(b) deriv2

(a) shaw

(b) deriv2

Figure 5.2: Relative errors of the regularized solutions computed by “s+d” and
“s+s” with respect to that by “d”, $\varepsilon=10^{-3}$.

First, we compare the relative errors $RE(k)$ for the three different
implementations of LSQR when $\varepsilon=10^{-3}$. From Fig. 5.1 we can find
the convergence behaviors of the three implementations “d”, “s+d” and “s+s”
are of highly consistence. The semi-convergence points $k_{0}$ are the same
and the relative error curves coincide until many steps after semi-
convergence, and thus the optimal regularized solutions computed by “d”, “s+d”
and “s+s” have the same accuracy. In order to give a more clear comparison
about accuracy of solutions, we also plot the relative error curves of $x_{k}$
computed by “s+d” and “s+s” with respect to that by “d”. Fig. 5.2 shows that
these two relative errors are much smaller than $RE(k)$ of “d” until semi-
convergence and this is also true for many iterations afterwards. These
results confirm that both the Lanczos bidiagonalization and updating procedure
can be implemented using single precision without sacrificing any accuracy of
regularized solutions for $\varepsilon=10^{-3}$.

(a) shaw

(b) deriv2

(c) gravity

(d) heat

Figure 5.3: Relative errors between regularized solutions computed by “s+d”
and “s+s”, $\varepsilon=10^{-3}$.

To further compare the accuracy of solutions computed by “s+d” and “s+s”, we
plot in Fig. 5.3 the relative error curves of these two solutions and their
upper bounds in (39). Here we set the upper bounds as
$\kappa(\widetilde{R}_{k})\mathbf{\bar{u}}$ with $\mathbf{\bar{u}}$ the
roundoff unit of single precision. From Fig. 5.3 we can find that
$\kappa(\widetilde{R}_{k})$ for the four test problems grow very slightly,
which lead to the upper bounds much smaller than $RE(k)$. Therefore, the
regularized solutions computed by “s+d” and “s+s” have the same accuracy,
which has already been clearly observed from Fig. 5.1.

Table 5.2: Comparison of relative errors $RE(k)$ and estimates of the optimal iterations $k_{0}$ by L-curve criterion and discrepancy principle ($\tau=1.001$), $\varepsilon=10^{-3}$. Work precision | shaw | deriv2 | gravity | heat
---|---|---|---|---
| Optimal |
d |  $0.0396$ ($8$) |  $0.1471$ ($15$) |  $0.0105$ ($10$) |  $0.0206$ ($22$)
s+d |  $0.0396$ ($8$) |  $0.1471$ ($15$) |  $0.0105$ ($10$) |  $0.0206$ ($22$)
s+s |  $0.0396$ ($8$) |  $0.1471$ ($15$) |  $0.0105$ ($10$) |  $0.0206$ ($22$)
| L-curve |
d |  $0.0396$ ($8$) |  $0.1529$ ($17$) |  $0.0105$ ($10$) |  $0.0392$ ($27$)
s+d |  $0.0396$ ($8$) |  $0.1529$ ($17$) |  $0.0105$ ($10$) |  $0.0392$ ($27$)
s+s |  $0.0396$ ($8$) |  $0.1529$ ($17$) |  $0.0105$ ($10$) |  $0.0392$ ($27$)
| Discrepancy principle |
d |  $0.0473$ ($7$) |  $0.1557$ ($13$) |  $0.0166$ ($8$) |  $0.0280$ ($19$)
s+d/s |  $0.0473$ ($7$) |  $0.1557$ ($13$) |  $0.0166$ ($8$) |  $0.0280$ ($19$)

Table 5.2 shows the relative errors of the regularized solutions at the semi-
convergence point $k_{0}$ and the estimates of $k_{0}$ by L-curve criterion
and discrepancy principle, where the corresponding iteration number is in
brackets. We find that the three optimal iterations $k_{0}$ and corresponding
$RE(k_{0})$ for “d”, “s+d” and “s+s” are the same, which has also been
observed from Fig. 5.1. This is also true for the L-curve criterion, and the
method gets an over-estimate of $k_{0}$ for heat. For the discrepancy
principle, by (41) we know that the estimates of $k_{0}$ for “s+d” and “s+s”
are always the same since they compute the same $\bar{\phi}_{k+1}$, and we
find these estimates are the same at that for “d”. The discrepancy principle
gets under-estimates of $k_{0}$ for the four test problems and thus the
solutions are over-regularized.

(a) shaw

(b) deriv2

(c) gravity

(d) heat

Figure 5.4: Semi-convergence curves for LSQR implemented using different
computing precisions, $\varepsilon=10^{-5}$.

(a) gravity

(b) heat

Figure 5.5: Semi-convergence curves for LSQR implemented using different
computing precisions, $\varepsilon=10^{-7}$.

From the above experimental results, we can make sure that single precision is
enough for LSQR for solving these four linear ill-posed problems when
$\varepsilon=10^{-3}$. We have also made numerical experiments for
$\varepsilon=10^{-4},10^{-5}$ and get similar results. Here we only show the
semi-convergence curves for $\varepsilon=10^{-5}$ in Fig. 5.4. For this noise
level, single precision is also enough for LSQR for solving the four test
problems. Note from subfigure (a) that the relative errors for “s+d” and “s+s”
quickly become bigger than that for “d” after semi-convergence. This reminds
us that for shaw, if $\varepsilon$ is smaller than $10^{-5}$, simply using a
single precision variant of LSQR will lead to a loss of accuracy of
regularized solutions. For extremely small noise level $\varepsilon=10^{-7}$,
we plot the semi-convergence curves for gravity and heat in Fig. 5.5. We can
clearly find that neither “s+d” nor “s+s” can compute a regularized solution
with accuracy as good as the best one achieved by “d”. In real applications,
the noise levels are seldom so small, and it is almost always possible to
implement LSQR using single precision to compute a regularized solution
without sacrificing accuracy.

###  5.2 Two dimensional case

We generate two image deblurring problems using codes from [11] for testing
two dimensional linear ill-posed problems, with the goal to restore an image
from a blurred and noisy one $b=Ax_{ex}+e$, where $x_{ex}$ denotes the true
image and $A$ denotes the blurring operator. For the background of image
deblurring, we refer the readers to [20]. For PRblurspeckle, which simulates
spatially invariant blurring caused by atmospheric turbulence, we use the true
image “Hubble Space Telescope” with image size of $N=128$ (i.e., the true and
blurred images have $128\times 128$ pixels), and the blur level is set to be
medium. For PRblurdefocus, which simulates a spatially invariant, out-of-focus
blur, we use the true image “Cameraman” with image size of $N=256$, and the
blur level is set to be severe. We use zero boundary condition for both the
two blurs to construct $A$. The two true images are shown in Fig. 5.6.

(a) Hubble Space Telescope

(b) Cameraman

Figure 5.6: True images for testing PRblurspeckle and PRblurdefocus.

For PRblurspeckle with image “Hubble Space Telescope”, the noise levels are
set as $\varepsilon=10^{-2}$ and $\varepsilon=10^{-3}$. Fig. 5.7 depicts the
semi-convergence curves for “d”, “s+d” and “s+s” as well as relative error
curves of $x_{k}$ computed by “s+d” and “s+s” with respect to that by “d”. For
both the two noise levels, we find that the curves of $RE(k)$ coincide until
many steps after semi-convergence and the three semi-convergence points are
the same. The relative errors of $x_{k}$ computed by “s+d” and “s+s” with
respect to that by “d” are shown in subfigures (c) and (d), and they are much
smaller than $RE(k)$ of “d” until many steps after semi-convergence. The
numerical results confirm that for PRblurspeckle, the LSQR can be implemented
using single precision without sacrificing accuracy of regularized solutions
for $\varepsilon=10^{-2}$ or $\varepsilon=10^{-3}$. Fig. 5.8 shows the blurred
images and corresponding restored ones for the two noise levels, where the
restored images are obtained from $x_{k}$ at the semi-convergence point (
$k_{0}$ and $RE(k_{0})$ are the same for “d”, “s+d” and “s+s”, and we choose
the $x_{k_{0}}$ computed by “s+s”). The result shows a good deblurring effect
of LSQR implemented using single precision. For noise levels smaller than
$10^{-4}$, we find that the semi-convergence behavior does not appear, which
means that noise amplification is tolerable even without regularization, and
thus we need not test for smaller noise cases.

(a) $\varepsilon=10^{-2}$

(b) $\varepsilon=10^{-2}$

(c) $\varepsilon=10^{-3}$

(d) $\varepsilon=10^{-3}$

Figure 5.7: Semi-convergence curves for LSQR implemented using different
computing precisions, and relative error curves of the regularized solutions
computed by “s+d” and “s+s” with respect to that by “d”, PRblurspeckle.

(a) $\varepsilon=10^{-2}$

(b) $\varepsilon=10^{-2}$

(c) $\varepsilon=10^{-3}$

(d) $\varepsilon=10^{-3}$

Figure 5.8: Images “Hubble Space Telescope” blurred by PRblurspeckle and
restored by LSQR at the optimal iteration: (a), (b) $\varepsilon=10^{-2}$;
(c), (d) $\varepsilon=10^{-3}$.

For PRblurdefocus with image “Cameraman”, the noise levels are set as
$\varepsilon=10^{-3}$ and $\varepsilon=10^{-4}$. For noise levels smaller than
$10^{-5}$, the semi-convergence behavior will not appear and we need not test
for those cases. Fig. 5.9 and Fig. 5.10 show the relative error curves and
blurred and restored images. Subfigures (b) and (d) of Fig. 5.9 depict the
relative errors between regularized solutions computed by “s+d” and “s+s” with
upper bounds $\kappa(\widetilde{R}_{k})\mathbf{\bar{u}}$. We can find that
$\kappa(\widetilde{R}_{k})$ for the two noise levels grow very slightly, which
lead to the upper bounds much smaller than $RE(k)$. Therefore, the regularized
solutions computed by “s+d” and “s+s” have the same accuracy, which can
clearly observed from subfigures (a) and (c). The other experimental results
are similar to those of PRblurdefocus and we do not illustrate them in detail
any longer.

(a) $\varepsilon=10^{-3}$

(b) $\varepsilon=10^{-3}$

(c) $\varepsilon=10^{-4}$

(d) $\varepsilon=10^{-4}$

Figure 5.9: Semi-convergence curves for LSQR implemented using different
computing precisions, and relative errors between regularized solutions
computed by “s+d” and “s+s”, PRblurdefocus.

(a) $\varepsilon=10^{-3}$

(b) $\varepsilon=10^{-3}$

(c) $\varepsilon=10^{-4}$

(d) $\varepsilon=10^{-4}$

Figure 5.10: Images “Cameraman” blurred by PRblurdefocus and restored by LSQR
at the optimal iteration: (a), (b) $\varepsilon=10^{-3}$; (c), (d)
$\varepsilon=10^{-4}$.

Table 5.3 shows the relative errors of the regularized solutions at the semi-
convergence point $k_{0}$ and the estimates of $k_{0}$ by L-curve criterion
and discrepancy principle. For PRblurspeckle, the three optimal iterations
$k_{0}$ and corresponding $RE(k_{0})$ for “d”, “s+d” and “s+s” are the same,
while for PRblurdefocus the $k_{0}$ for “d” and “s+d/s” are differed only by
one and the corresponding $RE(k_{0})$ are the same. For the discrepancy
principle, the estimates of $k_{0}$ are the same for “d” and “s+d/s”, and the
corresponding $RE(k_{0})$ is only slightly different for PRblurdefocus with
$\varepsilon=10^{-3}$. For the L-curve criterion, the estimates of $k_{0}$ and
the corresponding $RE(k_{0})$ for “d”, “s+d” and “s+s” are not the same for
some cases, but the differences are very slight.

From the above experimental results, we can make sure that single precision is
enough for LSQR for solving the 2-D image deblurring ill-posed problems if the
noise level is not extremely small. For those large scale problems, using
single precision instead of double precision can greatly improve the computing
efficiency while can still maintain the same accuracy for regularized
solutions.

Table 5.3: Comparison of relative errors $RE(k)$ and estimates of the optimal iterations $k_{0}$ by L-curve criterion and discrepancy principle ($\tau=1.001$). Work precision | PRblurspeckle | PRblurdefocus
---|---|---
| $\varepsilon=10^{-2}$ | $\varepsilon=10^{-3}$ | $\varepsilon=10^{-3}$ | $\varepsilon=10^{-4}$
| Optimal |
d |  $0.1850$ ($34$) |  $0.1102$ ($135$) |  $0.1058$ ($80$) |  $0.0542$ ($190$)
s+d |  $0.1850$ ($34$) |  $0.1102$ ($135$) |  $0.1058$ ($81$) |  $0.0542$ ($191$)
s+s |  $0.1850$ ($34$) |  $0.1102$ ($135$) |  $0.1058$ ($81$) |  $0.0542$ ($191$)
| L-curve |
d |  $0.1992$ ($47$) |  $0.1105$ ($149$) |  $0.1293$ ($96$) |  $0.0764$ ($100$)
s+d |  $0.1992$ ($47$) |  $0.1105$ ($149$) |  $0.1157$ ($91$) |  $0.0771$ ($100$)
s+s |  $0.1992$ ($47$) |  $0.1103$ ($145$) |  $0.1157$ ($91$) |  $0.0771$ ($100$)
| Discrepancy principle |
d |  $0.1937$ ($24$) |  $0.1200$ ($78$) |  $0.1081$ ($74$) |  $0.0599$ ($137$)
s+d/s |  $0.1937$ ($24$) |  $0.1200$ ($78$) |  $0.1099$ ($74$) |  $0.0599$ ($137$)

##  6 Conclusion

For the most commonly used iterative regularization algorithm LSQR for solving
linear discrete ill-posed problems, we have investigated how to exploit lower
precision computing in the two main parts of the algorithm, including the
construction of Lanczos vectors and updating procedure of iterative solutions.
Under an ideal model describing the ill-posed problem (1), we have shown that,
for not extremely small noise levels, single precision is enough for computing
Lanczos vectors with full reorthogonalization without loss of any accuracy of
the final regularized solution. For the procedure of updating $x_{k}$, we have
shown that the computation of $x_{i}$ and $w_{i}$, which is the most time
consuming part, can be performed using single precision without sacrificing
any accuracy as long as $\kappa(\widetilde{R}_{k})$ is not very big that is
almost always satisfied. We have made several numerical experiments to
illustrate numerical performances of the single precision variants of LSQR for
various types of discrete linear ill-posed problems.

Our results indicate that the most time consuming parts of the algorithm can
be implemented using single precision, and thus the performance of LSQR for
solving linear discrete ill-posed problems can be significantly enhanced. The
results can guide us towards an implementation that is both robust and
efficient. Future work includes developing practical C codes on high
performance computing architectures for specific real applications.

## Acknowledgments

The authors thank Dr. Long Wang and Professor Weile Jia for many useful
discussions about lower and mixed precision computations on high performance
computing architectures.

## References

  * [1] K. Ahmad, H. Sundar, and M. W. Hall, Data-driven mixed precision sparse matrix vector multiplication for GPUs, ACM Trans. Archit. Code Optim., 16 (2019), pp. 51:1–51:24, https://doi.org/10.1145/3371275.
  * [2] P. Amestoy, O. Boiteau, A. Buttari, M. Gerest, F. Jézéquel, J.-Y. l’Excellent, and T. Mary, Mixed precision low rank approximations and their application to block low rank LU factorization, (2021), https://hal.archives-ouvertes.fr/hal-03251738v1.
  * [3] Å. Björck, A bidiagonalization algorithm for solving large and sparse ill-posed systems of linear equations, BIT Numer. Math., 28 (1988), pp. 659–670.
  * [4] Å. Björck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia, 1996.
  * [5] P. Blanchard, N. J. Higham, F. Lopez, T. Mary, and S. Pranesh, Mixed precision block fused multiply-add: error analysis and application to GPU tensor cores, SIAM J. Sci. Comput., 42 (2020), pp. C124–C141, https://doi.org/10.1137/19M1289546.
  * [6] L. S. Borges, F. S. V. Bazán, and M. C. C. Cunha, Automatic stopping rule for iterative methods in discrete ill-posed problems, Comp. Appl. Math., 34 (2015), pp. 1175–1197.
  * [7] E. Carson, N. J. Higham, and S. Pranesh, Three-precision GMRES-based iterative refinement for least squares problems, SIAM J. Sci. Comput., 42 (2020), pp. A4063–A4083, https://doi.org/10.1137/20M1316822.
  * [8] J. Chung and S. Gazzola, Computational methods for large-scale inverse problems: a survey on hybrid projection methods, arXiv:2105.07221v2, (2021), https://doi.org/10.48550/arXiv.2105.07221v2.
  * [9] J. Chung, J. G. Nagy, and D. P. O’Leary, A weighted-GCV method for Lanczos-hybrid regularization., Electr. Trans. Numer. Anal., 28 (2008), pp. 149–167.
  * [10] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Kluwer Academic Publishers, 2000.
  * [11] S. Gazzola, P. C. Hansen, and J. G. Nagy, IR Tools: A MATLAB package of iterative regularization methods and large-scale test problems, Numer. Algor., 81 (2019), pp. 773–811.
  * [12] S. Gazzola, P. Novati, and M. R. Russo, On Krylov projection methods and Tikhonov regularization, Electr. Trans. Numer. Anal., 44 (2015), pp. 83–123.
  * [13] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 4th ed., 2013.
  * [14] G. H. Golub and H. G. Wahba, Generalized cross-validation as a method for choosing a good ridge parameter, Technometrics, 21 (1979), pp. 215–223.
  * [15] S. Gratton, E. Simon, D. Titley-Peloquin, and P. Toint, Exploiting variable precision in GMRES, arXiv:1907.10550v2, (2019), https://arxiv.org/abs/arXiv:1907.10550v2.
  * [16] P. C. Hansen, Analysis of discrete ill-posed problems by means of the l-curve, SIAM Review, 34 (1992), pp. 561–580.
  * [17] P. C. Hansen, Rank-deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion, SIAM, Philadelphia, 1998.
  * [18] P. C. Hansen, Regularization Tools version 4.0 for Matlab 7.3, Numer. Algor., 46 (2007), pp. 189–194.
  * [19] P. C. Hansen, Discrete Inverse Problems: Insight and Algorithms, SIAM, Philadelphia, 2010.
  * [20] P. C. Hansen, J. G. Nagy, and D. P. O’Leary, Deblurring Images: Matrices, Spectra and Filtering, SIAM, Philadelphia, 2006.
  * [21] N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, 2002.
  * [22] N. J. Higham, S. Pranesh, and M. Zounon, Squeezing a matrix into half precision, with an application to solving linear systems, SIAM J. Sci. Comput., 41 (2019), pp. A2536–A2551, https://doi.org/10.1137/18M1229511.
  * [23] I. Hnětynková, M. Plešinger, and Z. Strakoš, The regularizing effect of the Golub-Kahan iterative bidiagonalization and revealing the noise level in the data, BIT Numer. Math., 49 (2009), pp. 669–696, https://doi.org/10.1007/s10543-009-0239-7.
  * [24] J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, 2006.
  * [25] M. E. Kilmer and D. P. O’Leary, Choosing regularization parameters in iterative methods for ill-posed problems, SIAM J. Matrix Anal. Appl., 22 (2001), pp. 1204–1221.
  * [26] R. M. Larsen, Lanczos bidiagonalization with partial reorthogonalization, DAIMI Report Series, (1998).
  * [27] H. Li, G. Tan, and T. Zhao, Backward error analysis of the Lanczos bidiagonalization with reorthogonalization, arXiv:2210.10297v1, (2022), https://doi.org/10.48550/arXiv.2210.10297.
  * [28] F. Lopez and T. Mary, Mixed precision LU factorization on GPU tensor cores: Reducing data movement and memory footprint, (2020), http://eprints.maths.manchester.ac.uk/2782/.
  * [29] G. Meurant and Z. Strakoš, The Lanczos and conjugate gradient algorithms in finite precision arithmetic, Acta Numerica, 15 (2006), pp. 471–542.
  * [30] V. A. Morozov, Regularization of incorrectly posed problems and the choice of regularization parameter, USSR Computational Mathematics and Mathematical Physics, 6 (1966), pp. 242–251.
  * [31] F. Natterer, The Mathematics of Computerized Tomography, SIAM, Philadelphia, 2001.
  * [32] C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Software, 8 (1982), pp. 43–71.
  * [33] L. Reichel, H. Sadok, and W. H. Zhang, Simple stopping criteria for the LSQR method applied to discrete ill-posed problems, Numer. Algor., 84 (2020), pp. 1381–1395.
  * [34] R. A. Renaut, S. Vatankhah, and V. E. Ardesta, Hybrid and iteratively reweighted regularization by unbiased predictive risk and weighted GCV for projected systems, SIAM J. Sci. Statist. Comput., 39 (2017), pp. B221–B243.
  * [35] M. Richter, Inverse Problems: Basics, Theory and Applications in Geophysics, Springer, 2016, https://doi.org/10.1007/978-3-319-48384-9_3.
  * [36] A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems, Washington, DC, 1977.
  * [37] C. R. Vogel, Computational Methods for Inverse Problems, SIAM, Philadelphia, 2002.
  * [38] P.-Å. Wedin, Perturbation bounds in connection with singular value decomposition, BIT Numer. Math., 12 (1972), pp. 99–111.
  * [39] L. M. Yang, A. Fox, and G. Sanders, Rounding error analysis of mixed precision block householder QR algorithms, SIAM J. Sci. Comput., 43 (2021), pp. A1723–A1753, https://doi.org/10.1137/19M1296367.


|  $2,10,14,22,26,34$ |  $24\,(9Z_{A}-9Z_{D}-2Z_{new1}^{(34)}+2Z_{new2}^{(34)})$
|  $3,15,21,33$ |  $16\,(9Z_{A}+9Z_{D}+12Z_{new1}^{(34)}-2Z_{new2}^{(34)})$
|  $4,8,12,16,20,24,28,32$ |  $24\,(15Z_{A}-3Z_{D}-2Z_{new1}^{(34)}-2Z_{new2}^{(34)})$
|  $6,30$ |  $8\,(27Z_{A}-27Z_{D}+30Z_{new1}^{(34)}+2Z_{new2}^{(34)})$
|  $9,18,27$ |  $64\,(3Z_{new1}^{(34)}+Z_{new2}^{(34)})$
$(35,\frac{105}{37},148)$ |  $1,\dots,36$ |  $228\,Z_{A}$
$(36,\frac{54}{19},152)$ |  odd $,\neq 19$ |  $8\,(20Z_{A}+17Z_{D})$
|  even |  $160\,(2Z_{A}-Z_{D})$
|  $19$ |  $432\,Z_{D}$
$(37,\frac{37}{13},156)$ |  $\neq 0\mbox{ mod }3,13$ |  $12\,(14Z_{A}+5Z_{new}^{(37)})$
|  $=0\mbox{ mod }3$ |  $168\,(2Z_{A}-Z_{new}^{(37)})$
|  $=0\mbox{ mod }13$ |  $288\,Z_{new}^{(37)}$
$(38,\frac{57}{20},160)$ |  odd, $\neq 5,15,25,35$ |  $32\,(4Z_{A}+4Z_{D}+Z_{new1}^{(38)})$
|  $2,6,14,18,22,26,34,38$ |  $16\,(16Z_{A}-8Z_{D}+Z_{new1}^{(38)}+6Z_{new2}^{(38)})$
|  $4,12,28,36$ |  $64\,(8Z_{A}-4Z_{D}-Z_{new1}^{(38)}-3Z_{new2}^{(38)})$
|  $5,15,20,25,35$ |  $128\,(2Z_{A}+2Z_{D}-Z_{new1}^{(38)})$
|  $8,16,24,32$ |  $96\,(4Z_{A}-4Z_{D}+Z_{new1}^{(38)}-2Z_{new2}^{(38)})$
|  $10,30$ |  $64\,(2Z_{A}+2Z_{D}-Z_{new1}^{(38)}+6Z_{new2}^{(38)})$
$(39,\frac{117}{41},164)$ |  $1,\dots,40$ |  $252\,Z_{A}$
$(40,\frac{20}{7},168)$ |  odd, $\neq 3\mbox{ mod }6,\neq 7,35$ |  $8\,(16Z_{A}+40Z_{D}-15Z_{new1}^{(40)}+2Z_{new2}^{(40)})$
|  even, $\neq 0\mbox{ mod }6,\neq 14,28$ |  $32\,(8Z_{A}-4Z_{D}+Z_{new2}^{(40)})$
|  $3,9,15,27,33,39$ |  $16\,(16Z_{A}-8Z_{D}+9Z_{new1}^{(40)}-4Z_{new2}^{(40)})$
|  $6,12,18,24,30,36$ |  $128\,(4Z_{A}-2Z_{D}-Z_{new2}^{(40)})$
|  $7,35$ |  $48\,(3Z_{new1}^{(40)}+2Z_{new2}^{(40)})$
|  $14,28$ |  $192\,Z_{new2}^{(40)}$
|  $21$ |  $288\,Z_{new1}^{(40)}$
$(41,\frac{123}{43},172)$ |  $1,\dots,42$ |  $264\,Z_{A}$
|  |

In the above table, terms labelled $Z_{new}$ are unphysical modular
invariants. We list a few of them below. Notice that for $k=22$ there are two
unphysical invariants. This number grows for generic $k$.

$\begin{split}Z^{(13)}_{new}&=\absolutevalue{\chi_{1}-\chi_{11}}^{2}+\absolutevalue{\chi_{2}+\chi_{8}}^{2}+\absolutevalue{\chi_{4}-\chi_{14}}^{2}+\absolutevalue{\chi_{7}+\chi_{13}}^{2}+2\absolutevalue{\chi_{5}}^{2}+2\absolutevalue{\chi_{10}}^{2}\\\
\\\
Z^{(19)}_{new}&=\absolutevalue{\chi_{1}+\chi_{13}}^{2}+\absolutevalue{\chi_{2}-\chi_{16}}^{2}+\absolutevalue{\chi_{4}+\chi_{10}}^{2}+\absolutevalue{\chi_{5}-\chi_{19}}^{2}+\absolutevalue{\chi_{8}+\chi_{20}}^{2}+\absolutevalue{\chi_{11}+\chi_{17}}^{2}\\\
&\quad+2\absolutevalue{\chi_{7}}^{2}+2\absolutevalue{\chi_{14}}^{2}\\\ \\\
Z_{new}^{(18)}&=\absolutevalue{\chi_{1}-\chi_{9}}^{2}+\absolutevalue{\chi_{3}+\chi_{13}}^{2}+\absolutevalue{\chi_{4}+\chi_{16}}^{2}+\absolutevalue{\chi_{7}+\chi_{17}}^{2}+\absolutevalue{\chi_{8}+\chi_{12}}^{2}+\absolutevalue{\chi_{11}-\chi_{19}}^{2}\\\
\\\
Z_{new1}^{(22)}&=2\left(\absolutevalue{\chi_{1}-\chi_{17}}^{2}+\absolutevalue{\chi_{5}+\chi_{11}}^{2}+\absolutevalue{\chi_{7}-\chi_{23}}^{2}+\absolutevalue{\chi_{13}+\chi_{19}}^{2}+\absolutevalue{\chi_{8}+\chi_{16}}^{2}\right)\\\
&\quad+\absolutevalue{\chi_{2}+\chi_{10}+\chi_{14}+\chi_{22}}^{2}\\\
Z_{new2}^{(22)}&=4\absolutevalue{\chi_{6}-\chi_{18}}^{2}+2\absolutevalue{\chi_{2}+\chi_{10}-\chi_{14}-\chi_{22}}^{2}\\\
\end{split}$ (C.1)

##  Appendix D Sage calculations for unitary Virasoro minimal models

Table 8: Poincaré Sums for unitary Virasoro minimal models Minimal Model $(m,m^{\prime},c,N_{0})$ |  Seed Primary $(\chi_{r,s})$ |  Dimension $(h_{r,s})$ |  Partition function
---|---|---|---
Ising |  $\chi_{1,1},\,\chi_{1,2}$ |  $0,\frac{1}{2}$ |  $32\,Z_{AA}$
$(3,4,\frac{1}{2},16)$ |  $\chi_{2,2}$ |  $\frac{1}{16}$ |  $32\,Z_{AA}$
Tricritical Ising |  $\chi_{1,1},\,\chi_{1,3},\,\chi_{2,3},\,\chi_{3,3}$ |  $0,\frac{3}{2},\frac{3}{5},\frac{1}{10}$ |  $384\,Z_{AA}$
$(4,5,\frac{7}{10},80)$ |  $\chi_{1,2},\,\chi_{2,2}$ |  $\frac{7}{16},\frac{3}{80}$ |  $384\,Z_{AA}$
3-State Potts |  $\chi_{1,1},\,\chi_{1,2},\,\chi_{1,3},\,\chi_{1,4}$ |  $0,\frac{2}{5},\frac{7}{5},3$ |  $96(4Z_{AA}+Z_{AD})$
$(5,6,\frac{4}{5},120)$ |  $\chi_{3,3},\,\chi_{3,4}$ |  $\frac{1}{15},\frac{2}{3}$ |  $576\,Z_{AD}$
|  $\chi_{2,2},\,\chi_{2,3},\,\chi_{2,4},\,\chi_{4,4}$ |  $\frac{1}{40},\frac{21}{40},\frac{13}{8},\frac{1}{8}$ |  $384(2Z_{AA}-Z_{AD})$
Tricritical Potts |  $\chi_{1,1}$ |  $0$ |  $128(4Z_{AA}+Z_{AD})$
$(6,7,\frac{6}{7},168)$ |  $\chi_{1,3}$ |  $\frac{4}{3}$ |  $768\,Z_{AD}$
|  $\chi_{1,2}$ |  $\frac{3}{8}$ |  $512(2Z_{AA}-Z_{AD})$
Name |  $\chi_{1,1}$ |  $0$ |  $512(Z_{AA}+Z_{AD})$
$(7,8,\frac{25}{28},224)$ |  $\chi_{4,4}$ |  $\frac{15}{224}$ |  $512(Z_{AA}+Z_{AD})$
|  $\chi_{2,2}$ |  $\frac{3}{224}$ |  $256(7\,Z_{AA}-5\,Z_{AD})$
Name |  $\chi_{1,1}$ |  $0$ |  $572(Z_{AA}+Z_{AD})$
$(8,9,\frac{11}{12},288)$ |  $\chi_{1,4}$ |  $\frac{87}{32}$ |  $572(Z_{AA}+Z_{AD})$
|  $\chi_{1,2}$ |  $\frac{11}{32}$ |  $288(7\,Z_{AA}-5\,Z_{AD})$
|  $\chi_{3,3}$ |  $\frac{1}{36}$ |  $572(Z_{AA}+Z_{AD})$
|  $\chi_{3,4}$ |  $\frac{143}{288}$ |  $572(Z_{AA}+Z_{AD})$
|  $\chi_{3,6}$ |  $\frac{899}{288}$ |  $288(7\,Z_{AA}-5\,Z_{AD})$
Name |  $\chi_{1,1}$ |  $0$ |  $432(2Z_{AA}+Z_{AD})$
$(9,10,\frac{11}{12},360)$ |  $\chi_{2,2}$ |  $\frac{1}{120}$ |  $864(2Z_{AA}-Z_{AD})$
|  $\chi_{5,6}$ |  $\frac{28}{45}$ |  $1728Z_{AD}$
|  $\chi_{1,3}$ |  $\frac{11}{9}$ |  $432(2Z_{AA}+Z_{AD})$
|  $\chi_{2,3}$ |  $\frac{143}{360}$ |  $864(2Z_{AA}-Z_{AD})$
|  $\chi_{5,5}$ |  $\frac{1}{15}$ |  $1728\,Z_{AD}$
Name |  $\chi_{1,1}$ |  $0$ |  $576(2Z_{AA}+Z_{AD})$
$(10,11,\frac{52}{55},440)$ |  $\chi_{1,2}$ |  $\frac{13}{40}$ |  $1152(2Z_{AA}-Z_{AD})$
|  $\chi_{1,5}$ |  $\frac{23}{5}$ |  $2304\,Z_{AD}$
Name |  $\chi_{1,1}$ |  $0$ |  $768(Z_{AA}+Z_{AD}+Z^{(12)})$
$(11,12,\frac{21}{22},528)$ |  $\chi_{2,2}$ |  $\frac{1}{176}$ |  $768(4Z_{AA}-2Z_{AD}-Z^{(12)})$
|  $\chi_{3,3}$ |  $\frac{1}{66}$ |  $1536(Z_{AA}+Z_{AD}-Z^{(12)})$
|  $\chi_{4,4}$ |  $\frac{5}{176}$ |  $768(3Z_{AA}-3Z_{AD}+2Z^{(12)})$
|  $\chi_{6,6}$ |  $\frac{35}{528}$ |  $1536(Z_{AA}+Z_{AD}-Z^{(12)})$
|  |  |

## References

  * [1] A. Maloney and E. Witten, _Quantum Gravity Partition Functions in Three Dimensions_ , _JHEP_ 02 (2010) 029 [0712.0155].
  * [2] R. Dijkgraaf, J. M. Maldacena, G. W. Moore and E. P. Verlinde, _A Black Hole Farey Tail_ , hep-th/0005003.
  * [3] J. Manschot, _AdS(3) Partition Functions Reconstructed_ , _JHEP_ 10 (2007) 103 [0707.1159].
  * [4] J. Manschot and G. W. Moore, _A Modern Fareytail_ , _Communications in Number Theory and Physics_ 4 (2010) 103.
  * [5] S. Giombi, A. Maloney and X. Yin, _One-loop Partition Functions of 3D Gravity_ , _JHEP_ 08 (2008) 007 [0804.1773].
  * [6] C. A. Keller and A. Maloney, _Poincare Series, 3D Gravity and CFT Spectroscopy_ , _JHEP_ 02 (2015) 080 [1407.6008].
  * [7] A. Castro, M. R. Gaberdiel, T. Hartman, A. Maloney and R. Volpato, _The Gravity Dual of the Ising Model_ , _Phys. Rev._ D85 (2012) 024032 [1111.1987].
  * [8] C.-M. Jian, A. W. W. Ludwig, Z.-X. Luo, H.-Y. Sun and Z. Wang, _Establishing strongly-coupled 3D AdS quantum gravity with Ising dual using all-genus partition functions_ , _JHEP_ 10 (2020) 129 [1907.06656].
  * [9] P. Saad, S. H. Shenker and D. Stanford, _JT gravity as a matrix integral_ , 1903.11115.
  * [10] E. Witten, _Matrix Models and Deformations of JT Gravity_ , _Proc. Roy. Soc. Lond. A_ 476 (2020) 20200582 [2006.13414].
  * [11] A. Maloney and E. Witten, _Averaging over Narain moduli space_ , _JHEP_ 10 (2020) 187 [2006.04855].
  * [12] N. Afkhami-Jeddi, H. Cohn, T. Hartman and A. Tajdini, _Free partition functions and an averaged holographic duality_ , 2006.04839.
  * [13] P. Bantay, _The kernel of the modular representation and the Galois action in RCFT_ , _Commun. Math. Phys._ 233 (2003) 423 [math/0102149].
  * [14] P. Saad, S. H. Shenker and D. Stanford, _A semiclassical ramp in SYK and in gravity_ , 1806.06840.
  * [15] J. Cotler and K. Jensen, _AdS 3 gravity and random CFT_, 2006.08648.
  * [16] J. Cotler and K. Jensen, _AdS 3 wormholes from a modular bootstrap_, _JHEP_ 11 (2020) 058 [2007.15653].
  * [17] N. Benjamin, H. Ooguri, S.-H. Shao and Y. Wang, _Light-cone modular bootstrap and pure gravity_ , _Phys. Rev. D_ 100 (2019) 066029 [1906.04184].
  * [18] N. Benjamin, S. Collier and A. Maloney, _Pure Gravity and Conical Defects_ , _JHEP_ 09 (2020) 034 [2004.14428].
  * [19] L. F. Alday and J.-B. Bae, _Rademacher Expansions and the Spectrum of 2d CFT_ , _JHEP_ 11 (2020) 134 [2001.00022].
  * [20] H. Maxfield and G. J. Turiaci, _The path integral of 3D gravity near extremality; or, JT gravity with defects as a matrix integral_ , 2006.11317.
  * [21] A. Cappelli, C. Itzykson and J. B. Zuber, _Modular Invariant Partition Functions in Two-Dimensions_ , _Nucl. Phys._ B280 (1987) 445.
  * [22] A. Cappelli, C. Itzykson and J. B. Zuber, _The ADE Classification of Minimal and A1(1) Conformal Invariant Theories_ , _Commun. Math. Phys._ 113 (1987) 1.
  * [23] A. R. Chandra and S. Mukhi, _Towards a Classification of Two-Character Rational Conformal Field Theories_ , _JHEP_ 04 (2019) 153 [1810.09472].
  * [24] A. R. Chandra and S. Mukhi, _Curiosities above c = 24_ , _SciPost Phys._ 6 (2019) 053 [1812.05109].
  * [25] S. Mukhi, R. Poddar and P. Singh, _Rational CFT with three characters: the quasi-character approach_ , _JHEP_ 05 (2020) 003 [2002.01949].
  * [26] J. Bagger, N. Lambert, S. Mukhi and C. Papageorgakis, _Multiple Membranes in M-theory_ , _Phys. Rept._ 527 (2013) 1 [1203.3546].
  * [27] W. Stein et al., _Sage Mathematics Software (Version 2020.2)_.  The Sage Development Team, 2020.
  * [28] T. Gannon, _The Classification of affine SU(3) modular invariant partition functions_ , _Commun. Math. Phys._ 161 (1994) 233 [hep-th/9212060].
  * [29] T. Gannon, _The Classification of SU(3) modular invariants revisited_ , _Ann. Inst. H. Poincare Phys. Theor._ 65 (1996) 15 [hep-th/9404185].
  * [30] C. Itzykson, _Level One Kac-Moody Characters And Modular Invariance_ , in _Annecy 1988, Proceedings, Conformal Field Theories And Related Topics_ , 1988.
  * [31] P. Degiovanni, _Z/NZ Conformal Field Theories_ , _Commun. Math. Phys._ 127 (1990) 71.
  * [32] M. Bauer and C. Itzykson, _Modular transformations of SU(N) affine characters and their commutant_ , _Commun. Math. Phys._ 127 (1990) 617.
  * [33] A. N. Schellekens, _Meromorphic c = 24 conformal field theories_ , _Commun. Math. Phys._ 153 (1993) 159 [hep-th/9205072].
  * [34] P. Goddard and D. I. Olive, _Algebras, Lattices and Strings_ , in _Goddard, P. (Ed.), Olive, D. (Ed.): Ka-Moody And Virasoro Algebras* 210-255_ , 1983.
  * [35] A. Karch, Z.-X. Luo and H.-Y. Sun, _Holographic duality for Ising CFT with boundary_ , 2012.02067.
  * [36] L. Randall and R. Sundrum, _An Alternative to compactification_ , _Phys. Rev. Lett._ 83 (1999) 4690 [hep-th/9906064].
  * [37] J. L. Cardy, _Effect of Boundary Conditions on the Operator Content of Two-Dimensional Conformally Invariant Theories_ , _Nucl. Phys. B_ 275 (1986) 200.
  * [38] S. D. Mathur, S. Mukhi and A. Sen, _On the Classification of Rational Conformal Field Theories_ , _Phys. Lett._ B213 (1988) 303.
  * [39] S. Mukhi, _Classification of RCFT from Holomorphic Modular Bootstrap: A Status Report_ , in _Pollica Summer Workshop 2019: Mathematical and Geometric Tools for Conformal Field Theories (PSW2019) Pollica, Salerno, Italy, June 3-21, 2019_ , 2019, 1910.02973.