# Solving ill-conditioned linear algebraic systems using methods that improve
conditioning

A.S. Leonov1

1National Research Nuclear University MEPhI,
Moscow, Russian Federation
E-mail<EMAIL_ADDRESS>
> We consider the solution of systems of linear algebraic equations (SLAEs)
> with an ill-conditioned or degenerate exact matrix and an approximate right-
> hand side. An approach to solving such a problem is proposed and justified,
> which makes it possible to improve the conditionality of the SLAE matrix
> and, as a result, obtain an approximate solution that is stable to
> perturbations of the right hand side with higher accuracy than using other
> methods. The approach is implemented by an algorithm that uses so-called
> minimal pseudoinverse matrices. The results of numerical experiments are
> presented that confirm the theoretical provisions of the article.
>
> 2010 Mathematics subject classifications: 15A09, 15A10, 65F22, 15A29
>
> Keywords: ill-conditioned SLAEs, minimal pseudoinverse matrix.

##  1 Introduction

When studying various scientific problems, it is often necessary to solve
systems of linear algebraic equations (SLAEs) of very large dimensions. Such
systems often turn out to be ill-conditioned or degenerate. As a consequence,
classical methods of linear algebra for solving these SLAEs (see, for example,
[1, 2, 3, 4]) sometimes turn out to be of little use due to instability when
rounding errors accumulate. This problem can be partially solved by using
calculations on computers with double, quadruple and higher precision.
However, this approach is not applicable to all computing systems on personal
computers (PCs). Therefore, the practical solution of ill-conditioned and
degenerate SLAEs requires the use of special methods customized to combat
numerical instability.

Currently, many such methods are known with varying complexity and scope of
applicability (see [3, 5, 6], etc.). Let us especially mention here methods
that are relevant to the subsequent presentation. Among them, the
regularization method of A.N. Tikhonov and its variant in the form of the
regularized least squares method [7] are widely known. The TSVD (truncated
singular value decomposition) method in its various modifications is also
often used (see, for example, [8, 9, 10]). A minimal pseudoinverse matrix
(MPM) method has also been developed (see [11] and [12, 13]) that will be
mainly used further.

A distinctive feature of this work is the construction of special methods for
solving SLAEs, _significantly improving the condition number_ of their matrix
(i.e., significantly reducing them). This is how the proposed methods differ
from those listed above, in which reducing the condition number is not the
goal, but turns out to be a kind of by-product.

The article is organized as follows. In Sect.2, methods for stable solution of
SLAEs, which are used to varying degrees in the work, are briefly described.
In Sect.3, some auxiliary statements are considered. The problem statement,
the proposed solution scheme and its justification are presented in Sect.4. A
specific method that implements the proposed scheme for solving SLAEs is
proposed in Sect.5. The results of numerical experiments with model problems
and conclusions are presented in Sect.6.

##  2 On some methods for stable solution of SLAEs

Let $\mathbb{R}^{n}_{c},\mathbb{R}^{m}_{c}$ be spaces of column vectors with
Euclidean norms. We will solve SLAE of the form

$Az=u,\,\ z\in\mathbb{R}^{n}_{c},\,u\in\mathbb{R}^{m}_{c}$ (1)

with a real matrix $A,\,\|A\|\neq 0$, $\dim A=m\times n$, and right-hand side
$u\neq 0$. Since this system may not be solvable in in the classical sense, we
will consider its pseudosolutions $z^{\ast}\in\mathbb{R}^{n}_{c}$, i.e.
elements satisfying equality

$\left\|Az^{\ast}-u\right\|=\inf\\{\left\|Az-u\right\|:\,\
z\in\mathbb{R}^{n}_{c}\\}.$

We denote the set of all pseudo-solutions for the exact data of the problem
$\\{A,u\\}$ as $Z^{\ast}$. For a compatible system, these pseudo-solutions
represent conventional solutions. For any system of the form (1) there is a
unique normal pseudosolution $\bar{z}$, i.e. such a pseudo-solution for which

$\left\|\bar{z}\right\|=\inf\\{\left\|z\right\|:\,z\in Z^{\ast}\\}.$

For a solvable system (1) it is called a normal solution. For a uniquely
solvable system, the normal solution coincides with the classical solution.
Our goal is to find the vector $\bar{z}$ from the data $\\{A,u\\}$. The formal
solution to this problem is given using the pseudo-inverse matrix $A^{+}$:
$\bar{z}={A^{+}}u$ (see, for example, [14, 3, 4]).

Often this formulation of the problem is modified. It is believed that instead
of exact data $\\{A,u\\}$ we have at our disposal their approximations
$\\{A_{h},u_{\delta}\\}$, $\mathrm{dim}A_{h}=\mathrm{dim}A$,
$u_{\delta}\in\mathbb{R}^{m}_{c}$, and error levels, i.e. numbers
$\eta=(h,\delta)$, such that $\left\|A_{h}-A\right\|\leq
h,\,\left\|u_{\delta}-u\right\|\leq\delta$. Here and in what follows, the
norms of matrices are considered Euclidean:
$\left\|A\right\|={\left\|A\right\|_{E}}$. It is required to construct
approximations $z_{\eta}\in\mathbb{R}^{n}_{c}$ to the solution $\bar{z}$,
which have the stability property: $\left\|z_{\eta}-\bar{z}\right\|\rightarrow
0$ at $\eta\rightarrow 0$.

The solution to this problem using the formula
${z_{\eta}}=A_{h}^{+}{u_{\delta}}$ turns out to be unstable in the general
case, because convergence of $A_{h}^{+}\to{A^{+}}$ for $h\to 0$ is not
guaranteed. Therefore, other methods have been proposed. For example, we can
take as a solution the element $z_{\eta}=z_{\eta}^{\alpha(\eta)}$, which is
obtained by A.N. Tikhonovs regularization method:

$z_{\eta}^{\alpha(\eta)}={(\alpha(\eta)I+A_{h}^{*}{A_{h}})^{-1}}A_{h}^{*}{u_{\delta}}=T_{1}(A_{h})u_{\delta}.$

Here the parameter $\alpha(\eta)>0$ is chosen by algorithms described in [7,
12, 13, 15, 16] and other works. These algorithms provide the required
convergence of approximations. To implement the regularization method, it is
necessary to know the errors levels of the matrix and the right-hand side of
the SLAE.

Another solution to the problem can be obtained using the TSVD method. It is
based on the artificial zeroing of small singular values of the system matrix.
The method implementation scheme is as follows. The singular value
decomposition (SVD) of the approximate matrix is found:

$A_{h}=U_{h}R_{h}V^{T}_{h},\,R_{h}=\mathrm{diag}(\rho_{1},\rho_{2},...,\rho_{M}),\,\rho_{1}\geq\rho_{2}\geq...\geq\rho_{M}\geq
0$

with $M=\min(m,n)$. Next, using some method, the number $\kappa\in\mathbb{N}$
is selected that plays the role of a regularization parameter for which
$\rho_{\kappa}>0$ and $\rho_{k}=0$ for $k=\kappa+1,\kappa+2,...,M$. Then a new
matrix of singular numbers
$R_{\kappa}=\mathrm{diag}(\rho_{1},...,\rho_{\kappa},0,...,0)$ is constructed.
From it, as well as from the orthogonal matrices $U_{h},V_{h}$ of the SVD, an
approximate solution of the SLAE is found:

$z_{\kappa}=V_{h}R_{\kappa}^{+}U^{T}_{h}u_{\delta}=T_{2}(A_{h})u_{\delta},\,\
R_{\kappa}^{+}=\mathrm{diag}(\rho_{1}^{-1},...,\rho_{\kappa}^{-1},0,...,0).$

This corresponds to using, instead of the matrix $A_{h}$, another approximate
matrix $\hat{A}_{h}=U_{h}R_{\kappa}V^{T}_{h}$ and its pseudo-inversion. The
question of what is considered a small singular number, i.e. how to choose
$\kappa$ is decided differently. For example, in the article [10] it is
proposed to choose $\kappa=\kappa(h)$ so that

$(\rho_{\kappa+1}^{2}+\rho_{\kappa+2}^{2}+...+\rho_{M}^{2})^{1/2}\leq h,\,\
(\rho_{\kappa}^{2}+\rho_{\kappa+1}^{2}+...+\rho_{M}^{2})^{1/2}>h.$

This corresponds to choosing from all matrices $R_{\kappa}$ satisfying the
condition $\left\|R_{h}-R_{\kappa}\right\|_{E}\leq h$ the matrix that has the
minimum rank. This method requires knowledge of the matrix error $h$, but does
not use the $\delta$ value.

Both considered methods are based on stable approximation of the pseudoinverse
matrix $A^{+}$ by matrices $T_{1,2}(A_{h})$. Note that the regularization
method uses a full-rank matrix $T_{1}(A_{h})$, while the matrix $T_{2}(A_{h})$
of the TSVD method may have a significantly lower rank than the original
matrix $A_{h}$. The condition number of matrices $T_{1,2}(A_{h})$ is also
important. In what follows we use the spectral condition numbers
$\nu_{s}(A)=\rho_{1}/\rho_{r}$, where $r$ is the rank of the corresponding
matrix ${A}$. We will compare the values $\nu_{s}(T_{1,2}(A_{h}))$ and the
condition numbers of matrices arising in the methods proposed below.

In the works [11, 12, 13], a method of minimal pseudoinverse matrix (MPM) was
proposed and studied, which not only allows one to stably calculate $A^{+}$,
but also has a number of optimal properties. Let us briefly describe this
method. We introduce the space $\mathcal{A}$, consisting of matrices of size
$m\times n$ with Euclidean norms. Let $\bar{A}\in\mathcal{A}$ be an unknown
exact matrix for which we need to find a pseudo-inverse. Approximate data for
this are known quantities $A_{h}\in\mathcal{A}$ and $h$ such that
$\left\|A_{h}-\bar{A}\right\|\leq h$. Next, we define a set of matrices
$\mathcal{A}_{h}=\\{A\in\mathcal{A}:\,\left\|A_{h}-A\right\|\leq h\\}$,
comparable to the matrix $A_{h}$ in accuracy. It is clear that
$\bar{A}\in\mathcal{A}_{h}$. Consider the following extremal problem: find a
matrix $\tilde{A}_{h}\in\mathcal{A}$ for which equality

$\left\|\tilde{A}_{h}^{+}\right\|=\inf\\{\left\|A^{+}\right\|:\,A\in\mathcal{A}_{h}\\}$
(2)

is satisfied. The solution to the problem (2) is called the _matrix of the MPM
method_ , and the corresponding pseudoinverse matrix $\tilde{A}_{h}^{+}$ is a
_minimal pseudoinverse matrix_. It is proved that these matrices exist, and
for $h\rightarrow 0$ the accuracy estimate is valid:
$\left\|\tilde{A}_{h}^{+}-\bar{A}^{+}\right\|\thicksim
2h\left\|\bar{A}^{+}\right\|^{2}$. The estimate is asymptotically
unimprovable. This means that the MPM method gives the optimal order of
accuracy of approximation to the pseudoinverse matrix $\bar{A}^{+}$.

The MPM method is easily implemented numerically. The corresponding algorithm
consists of the following steps (see [11], [12, Chapter 5]).

1\. Finding the SVD of given matrix $A_{h}$: $A_{h}=U_{h}R_{h}V_{h}^{T}$,
where $R_{h}=\mathrm{diag}(\rho_{1},\rho_{2},...,\rho_{M})$, $\rho_{1}\geq\
rho_{2}\geq...\geq\rho_{M}\geq 0$, $M=\min(m,n)$, and $U_{h},V_{h}$ are
orthogonal matrices of sizes $m\times m$ and $n\times n$.

2\. A) Finding the quantities $\lambda_{k}=\frac{27}{16}\rho_{k}^{4}$ and
calculating the functions

$\rho_{k}(\lambda)=\\{\rho_{k}x_{k}(\lambda),\,0<\lambda\leq\lambda_{k};\,0,\,\lambda\geq\lambda_{k}\\};\,\rho_{k}(0)=\rho_{k},\,\,k=1,2,...,M,$

for $\lambda\geq 0$. Here $x_{k}(\lambda)$ is a solution to the equation
$x^{4}-x^{3}=\lambda\rho_{k}^{-4}$, belonging to the segment
$[1,\frac{3}{2}]$.

B) calculation of the function

$\beta(\lambda)=\sum_{k=1}^{M}[\rho_{k}(\lambda)-\rho_{k}]^{2},\,\lambda\geq
0.$

3\. Calculation of the solution $\lambda(h)>0$ of the equation
$\beta(\lambda)=h^{2}$. In practice, steps 2 A), B) are used repeatedly here.

4\. Finding a regularized singular value matrix

$\tilde{R}_{h}=\mathrm{diag}(\rho_{1}[\lambda(h)],\rho_{2}[\lambda(h)],...,\rho_{M}[\lambda(h)])$

and its pseudo-inverse

$\tilde{R}_{h}^{+}=\mathrm{diag}\\{\theta[\rho_{1}[\lambda(h)]],\theta[\rho_{2}[\lambda(h)]],...,\theta[\rho_{M}[\lambda(h)]]\\},$

where $\theta(\rho)=\\{\rho^{-1},\,\rho>0;\,0,\,\rho=0\\}$.

5\. The final calculation of the minimum pseudoinverse matrix:
$\tilde{A}_{h}^{+}=V_{h}\tilde{R}_{h}^{+}U_{h}^{T}$ and, if necessary, finding
the matrix of the MPM method $\tilde{A}_{h}=$ $U_{h}\tilde{R}_{h}V_{h}^{T}$.

As a commentary on the algorithm, we note that the function $\beta(\lambda)$
increases monotonically, having as its breakpoints the numbers $\lambda_{k}$.
At these points the uniqueness of the function $\beta(\lambda)$ is violated,
and at them it takes the values $\beta(\lambda_{k}-0)$ and
$\beta(\lambda_{k}+0)$. Therefore, when speaking about solving the equation
$\beta(\lambda)=h^{2}$, we mean finding a generalized solution, i.e. or jump
points of the function $\beta(\lambda)$ through the value $h^{2}$, or an
ordinary root. The issues of implementing the MPM algorithm, as well as the
influence of errors in calculating singular numbers for the original matrix on
the algorithm, are discussed in more detail in [11], [12]. It should also be
noted that the MPM method produces a matrix $\tilde{A}_{h}$ that is optimal in
order of the condition number.The method does not require knowledge of
$\delta$. Below we will use minimal pseudoinverse matrices in a new algorithm
for solving SLAEs, which guarantees an improvement in the condition number of
the system matrix.

##  3 Auxiliary statements

In this section we assume that for the exact matrix $\bar{A}\in\mathcal{A}$,
$\|\bar{A}\|\neq 0$, a family of approximate matrices
$\tilde{A}_{h}\in\mathcal{A}$, $\forall\,h\geq 0$ is given, such that the
approximation condition $\left\|\tilde{A}_{h}-\bar{A}\right\|\leq h$ is
satisfied.

###### Lemma 1

Let the matrices $\tilde{A}_{h}$ satisfy the relation
$\left\|\tilde{A}_{h}^{+}\right\|\leq\left\|\bar{A}^{+}\right\|$, $\forall
h,\leavevmode\nobreak\ 0\leq h<\left\|\bar{A}^{+}\right\|^{-1}$. Then: A)
$\mathrm{Rg}\tilde{A}_{h}=\mathrm{Rg}\bar{A}$; B)
$\left\|\tilde{A}_{h}^{+}-\bar{A}^{+}\right\|\leq
h\left\|\bar{A}^{+}\right\|^{2}(1-h\left\|\bar{A}^{+}\right\|)^{-3}$.

Proof. We denote as $\bar{\rho}_{k},\tilde{\rho}_{k},\leavevmode\nobreak\
1\leq k\leq M$ the singular values of the matrices $\bar{A},\tilde{A}_{h}$.
Then (see [11]) the relations

$\sum\limits_{k=1}^{M}(\tilde{\rho}_{k}-\bar{\rho}_{k})^{2}\leq\left\|\tilde{A}_{h}-\bar{A}\right\|^{2}\leq
h^{2},\leavevmode\nobreak\
\left\|\tilde{A}_{h}^{+}\right\|^{2}=\sum\limits_{k=1}^{{r}_{h}}\tilde{\rho}_{k}^{-2},\leavevmode\nobreak\
\left\|\bar{A}^{+}\right\|^{2}=\sum\limits_{k=1}^{r}\bar{\rho}_{k}^{-2},$ (3)

are satisfied, where $r_{h}=\mathrm{Rg}\tilde{A}_{h},r=\mathrm{Rg}\bar{A}$.
Let us prove A), i.e. equality $r_{h}=r$ for $0\leq
h<\left\|\bar{A}^{+}\right\|^{-1}$. We first assume that $r_{h}<r$. Then
$\tilde{\rho}_{k}=0$ for $k>r$, and therefore, taking into account (3) we
obtain

$\bar{\rho}_{r}\leq\sum\limits_{k=1}^{r-1}(\tilde{\rho}_{k}-\bar{\rho}_{k})^{2}+\sum\limits_{k=r}^{M}\bar{\rho}_{k}^{2}=\sum\limits_{k=1}^{M}(\tilde{\rho}_{k}-\bar{\rho}_{k})^{2}\leq
h^{2},\leavevmode\nobreak\
\bar{\rho}_{r}^{-2}\leq\left\|\bar{A}^{+}\right\|^{2}.$

It follows that $h^{2}\geq\left\|\bar{A}^{+}\right\|^{2}$, and this
contradicts the condition on $h$.

If we assume that $r_{h}>r$, then $\bar{\rho}_{r_{h}}=0$, and similarly we
obtain

$\tilde{\rho}_{r_{h}}^{2}=(\tilde{\rho}_{r_{h}}-\bar{\rho}_{r_{h}})^{2}\leq\sum\limits_{k=1}^{M}(\tilde{\rho}_{k}-\bar{\rho}_{k})^{2}\leq
h^{2},\leavevmode\nobreak\
\tilde{\rho}_{r_{h}}^{-2}\leq\sum\limits_{k=1}^{r_{h}}\tilde{\rho}_{k}^{-2}=\left\|\tilde{A}_{h}^{+}\right\|^{2}\leq\left\|\bar{A}^{+}\right\|^{2}.$

This again implies the inequality $h^{2}\geq\left\|\bar{A}^{+}\right\|^{2}$,
which contradicts the conditions. As a result we get that
$r=\mathrm{Rg}\bar{A}=r_{h}=\mathrm{Rg}\tilde{A}_{h}$.

The proof of part B) is based on the application of the estimate established
in [11]:

$\left\|A^{+}-B^{+}\right\|\leq\left\|A-B\right\|\left\|A^{+}\right\|^{2}/\left(1-\left\|A-B\right\|\left\|A^{+}\right\|\right)^{3}$

for matrices with $\dim A=\dim B,\leavevmode\nobreak\
\mathrm{Rg}A=\mathrm{Rg}B$. Setting $A=\bar{A},\leavevmode\nobreak\
B=\tilde{A}_{h}$ in this estimate and taking into account the proven part A)
of the lemma, we obtain the required inequality. $\square$

###### Corollary 1

If $\tilde{A}_{h}\rightarrow\bar{A}$ and
$\left\|\tilde{A}_{h}^{+}\right\|\leq\left\|\bar{A}^{+}\right\|$ for
$h\rightarrow 0$, then $\tilde{A}_{h}^{+}\rightarrow\bar{A}^{+}$.

This leads to next corollary.

###### Corollary 2

Let $\delta\to 0$, and for some function $h(\delta)\to 0$ a matrix
$A_{h(\delta)}$ is found such that

$\left\|A_{h(\delta)}-A\right\|\rightarrow
0,\,\left\|{A_{h(\delta)}^{+}}\right\|\leq\left\|{{A^{+}}}\right\|.$

Then $A_{h(\delta)}^{+}\to{A^{+}}$ for $\delta\to 0$ and therefore
${z_{\delta}}=A_{h(\delta)}^{+}{u_{\delta}}\to\bar{z}$.

##  4 Detailed formulation of the problem, proposed solution scheme and its
justification

In some practical cases, the SLAE matrix is known exactly, often in the form
of an analytical expression. In addition, it can be ill-conditioned or even
degenerate. Then there are reasons to replace it with a matrix close to it
with a better condition number and use the latter instead of the exact one.
This can be illustrated by the following example.

###### Example 1

Application of TSVD method to improve matrix conditioning.

In problems related to potential fields continuation, SLAE with a matrix of
the form

$\bar{A}=\left[\frac{1}{(x_{i}-y_{j})^{2}+H^{2}_{0}}\right]_{i,j=1}^{m,n},$
(4)

must sometimes be solved. Here $\\{x_{i}\\},\\{y_{j}\\}$ are uniform grids on
the segment $[-1,1]$. For $m=1991,n=2001$ and $H_{0}=0.1$ the matrix has a
condition number $\nu_{s}=\rho_{1}/\rho_{M}\approx 3.37\cdot 10^{19}$, i.e. is
extremely poorly conditioned. This is due to the specific (exponential) order
of decrease of the singular values of the matrix (see Fig.1A).

Figure 1: A) To Sect.4. Matrix singular values (4). B) To Sect.5. Solution
$h(\delta)$ of the equation
$\bar{\beta}_{\delta}(h)=\frac{b_{\delta}}{\left\|\bar{u}\right\|_{2}}$ c
$b_{\delta}=\sqrt{\mu_{\delta}^{2}+\delta^{2}}$ for SLAE with matrix (4).

By applying the TSVD method to the matrix, i.e., for example, replacing its
singular values that satisfy the condition $\rho_{k}<h$ with zeros, it is
possible, with an adequate choice of the value $h$, to improve the condition
number. Thus, for $h=10^{-10}$ we obtain a better condition number:
$\nu_{s}\approx 2.48\cdot 10^{9}$. Accordingly, the stability of numerical
solutions of SLAEs with such a matrix improves. However, it is not clear how
to constructively find $h$, since the proximity estimate for the matrices,
$h$, is unknown in the formulation under consideration. $\square$

We propose the following scheme for solving SLAEs using an approximate matrix
with an improved condition number. Let the SVD of the exact matrix be given:
$\bar{A}=U\bar{R}V^{T}$. Here $U,V$ are orthogonal matrices of size $m\times
m$ and $n\times n$, respectively, and $\bar{R}$ is a diagonal matrix of size
$m\times n$ containing the singular values of the matrix $\bar{A}$, sorted in
non-increasing order:

$\bar{R}=\mathrm{diag}\left[\bar{\rho}_{1},\bar{\rho}_{2},...,\bar{\rho}_{\bar{r}},0,...,0\right];\bar{\rho}_{k}\geq\bar{\rho}_{k+1}\geq
0,k=1,2,...,M-1;\,\bar{r}=\mathrm{Rg}\leavevmode\nobreak\ \bar{A}\leq M.$

We assume that the exact right-hand side $\bar{u}$ of the SLAE and the matrix
$\bar{A}$ are such that $\|\bar{u}\|\neq 0$, $\|\bar{A}\|\neq 0$. Then the
approximate right-hand side of the SLAE, $u_{\delta}$, satisfies the condition
$\|u_{\delta}\|\neq 0$ for sufficiently small values of $\delta$. It is these
$\delta$ that we consider further.

Scheme of solving the SLAE with data $\left\\{\bar{A},u_{\delta}\right\\}$.

1) For each $h,\leavevmode\nobreak\ 0\leq h\leq H=\mathrm{const}$, we look for
approximate matrices of the form $\tilde{A}_{h}=U\tilde{R}_{h}V^{T}$, where

$\tilde{R}_{h}=\mathrm{diag}\left[\tilde{\rho}_{1}x_{1}(h),\tilde{\rho}_{2}x_{2}(h),...,\tilde{\rho}_{\bar{r}}x_{\bar{r}}(h),0,...,0\right].$

The choice of functions $x_{k}(h),\,k=1,...,\bar{r}$ will be discussed below.
The number $H$ must be specified in each method under consideration.

2) We introduce the function
$\beta_{\delta}(h)=\left\|{\bar{A}\tilde{A}_{h}^{+}{u_{\delta}}-{u_{\delta}}}\right\|$,
$0\leq h\leq H$, and solve the equation
$\beta^{2}_{\delta}(h)=\delta^{2}+\beta^{2}_{\delta}(0)$. We denote its
solution as $h(\delta)>0$. Questions about the solvability of the equation
will be studied further.

3) We find the matrix $\tilde{A}_{h(\delta)}$ and use it to calculate the
approximate solution of the SLAE:
$z_{\delta}=\tilde{A}_{h(\delta)}^{+}u_{\delta}$.

Let’s make assumptions regarding the functions $x_{k}(h)$: for all $k=1,...,M$
the following requirements must be met.

A) $1<x_{k}(h)\leq c_{k}=\mathrm{const}$ $0<h\leq H$;

B) $x_{k}(+0)=x_{k}(0)=1,\leavevmode\nobreak\ x_{k}(H)=0$;

C) functions $x_{k}(h)$ are left continuous for all $h\in(0,H]$;

D) the functions $\theta[x_{k}(h)]$ are non-increasing as $0\leq h\leq H$.

Let us note some simple consequences from these assumptions.

E1) $x_{k}(h)=1\leavevmode\nobreak\ \Longleftrightarrow\leavevmode\nobreak\
h=0$; E2) $0\leq\theta[x_{k}(h)]\leq 1$; E3)
$\underset{h\rightarrow+0}{\lim}\theta[x_{k}(h)]=\theta[x_{k}(0)]=1$; E4) for
any $h_{0},\leavevmode\nobreak\ h_{0}\in[0,H]$, the equality holds

$\overline{\underset{h\rightarrow h_{0}}{\lim}}\leavevmode\nobreak\
\theta[x_{k}(h)]=\theta[x_{k}(h_{0})].$

We also note the inequality following from E2) and the definition of the
function $\theta$

$\left\|\tilde{A}_{h}^{+}\right\|^{2}=\sum\limits_{k=1}^{\bar{r}}\bar{\rho}_{k}^{-2}\theta[x_{k}^{2}(h)]\leq\sum\limits_{k=1}^{\bar{r}}\bar{\rho}_{k}^{-2}=\left\|\bar{A}^{+}\right\|^{2}.$
(5)

###### Theorem 1

Let conditions A) - D) be satisfied and, in addition, for each
$\delta,\leavevmode\nobreak\ 0<\delta<\delta_{0}=\mathrm{const}$, the
inequality
$\left\|u_{\delta}\right\|>\mu_{\delta}=\left\|\bar{A}\bar{A}^{+}u_{\delta}-u_{\delta}\right\|$
is valid. Then:

1) the function $\beta_{\delta}(h)$ non-decreasing as $h\in[0,H]$ and is left
continuous at every point $h>0$;

2) $\beta^{2}_{\delta}(+0)=\mu_{\delta}^{2},\leavevmode\nobreak\
\beta^{2}_{\delta}(\left\|\bar{A}\right\|)=\left\|u_{\delta}\right\|^{2}$;

3) the equation $\beta^{2}_{\delta}(h)=\delta^{2}+\mu^{2}_{\delta}$ has a
generalized solution $h(\delta)>0$;

4) $h(\delta)\rightarrow 0$ $\delta\rightarrow 0$.

###### Remark 1

In item 3) the equation with the monotonic function $\beta_{\delta}(h)$ is
solved. Its generalized solution $h(\delta)$ is a point for which the
inequalities hold

$\underset{h\rightarrow
h(\delta)}{\underline{\lim}}\beta^{2}_{\delta}(h)=\beta^{2}_{\delta}(h(\delta)-0)\leq\mu_{\delta}^{2}+\delta^{2}\leavevmode\nobreak\
\leq\beta^{2}_{\delta}(h(\delta)+0)=\underset{h\rightarrow
h(\delta)}{\overline{\lim}}\ \beta^{2}_{\delta}(h).$

Such equations arose in the study of the discrepancy principle for solving
ill-posed inverse problems (see, for example, [12]). In what follows we will
use the part of this inequality that takes into account the left continuity of
the function $\beta_{\delta}(h)$:

$\beta^{2}_{\delta}(h(\delta))=\beta^{2}_{\delta}(h(\delta)-0)\leq\delta^{2}+\mu^{2}_{\delta}.$
(6)

Proof. First, we prove 1). Equality

$\beta^{2}_{\delta}(h)=\left\|\bar{A}^{+}\tilde{A}_{h}^{+}u_{\delta}-u_{\delta}\right\|^{2}=\left\|U\bar{R}V^{T}V\tilde{R}_{h}^{+}U^{T}u_{\delta}-UU^{T}u_{\delta}\right\|^{2}=\left\|\bar{R}\tilde{R}_{h}^{+}v_{\delta}-v_{\delta}\right\|^{2},$

$v_{\delta}=U^{T}u_{\delta}=(v_{k})_{k=1}^{m}$ fulfilled. Next, taking into
account the form of the matrices $\bar{R}$ and $\tilde{R}_{h}^{+}$, we obtain

$\displaystyle\beta^{2}_{\delta}(h)=\sum\limits_{k=1}^{m}(\bar{\rho}_{k}\theta[\bar{\rho}_{k}x_{k}(h)]-1)^{2}v_{k}^{2}=\sum\limits_{k=1}^{\bar{r}}(\bar{\rho}_{k}\bar{\rho}_{k}^{-1}\theta[x_{k}(h)]-1)^{2}v_{k}^{2}+\sum\limits_{k=\bar{r}+1}^{m}v_{k}^{2}=$
$\displaystyle=\sum\limits_{k=1}^{\bar{r}}(1-\theta[x_{k}(h)])^{2}v_{k}^{2}+\mu_{\delta}^{2}.$
(7)

From this, as well as from properties D) and E2), the monotonicity (non-
decreasing) of the function $\beta_{\delta}(h)$ follows.

The equalities indicated in part 2) of the theorem follow from properties B),
E3) and E4). The existence of a generalized solution to the equation from item
3) follows from the monotonicity of the function $\beta_{\delta}(h)$, the
proven relations of item 2) and the condition
$\left\|u_{\delta}\right\|>\mu_{\delta}$ of theorem (for a similar statement,
see [11]).

Let us prove the convergence from part 4). Suppose that $h(\delta)\nrightarrow
0$ for $\delta\rightarrow 0$. Then there is a sequence
$\\{\delta_{N}\\},\delta_{N}\rightarrow 0$ for which
$h(\delta_{N})=h_{N}\rightarrow h_{0}>0$ as $N\rightarrow\infty$. From the
equality (7) and the property (6) of solutions $h(\delta_{N})=h_{N}$ it
follows:

$\beta^{2}_{\delta}(h_{N})=\sum\limits_{k=1}^{\bar{r}}(\theta[x_{k}(h_{N})]-1)^{2}(v_{\delta_{N}})_{k}^{2}+\mu_{\delta_{N}}^{2}\leq\delta_{N}^{2}+\mu_{\delta_{N}}^{2}\Longrightarrow\sum\limits_{k=1}^{\bar{r}}(\theta[x_{k}(h_{N})]-1)^{2}(v_{\delta_{N}})_{k}^{2}\leq\delta_{N}^{2}.$

Passing to the limit at $N\rightarrow\infty$ in the last inequality and taking
into account E4), we obtain:

$\sum\limits_{k=1}^{\bar{r}}(\theta[x_{k}(h_{0})]-1)^{2}\bar{v}_{k}^{2}=0.$

where $(\bar{v}_{k})_{k=1}^{m}=U^{T}\bar{u}$. It follows that there is a
number $k$ for which $\theta[x_{k}(h_{0})]=1$, because
$\left\|\bar{u}\right\|=\left\|U^{T}\bar{u}\right\|\neq 0$. But then, by
properties D), E1), E2) it turns out that $x_{k}(h_{0})=1\Longrightarrow
h_{0}=0$, i.e. $h_{N}\rightarrow 0$. The resulting contradiction proves the
theorem.

###### Theorem 2

The approximate solution of the SLAE,
$z_{\delta}=\tilde{A}_{h(\delta)}^{+}u_{\delta}$, converges to the normal
pseudo-solution $\bar{z}$ for $\delta\rightarrow 0$ .

Proof. From the form of the matrices $\tilde{A}_{h(\delta)},\bar{A}$ it
follows that

$\left\|\tilde{A}_{h(\delta)}-\bar{A}\right\|^{2}=\sum\limits_{k=1}^{\bar{r}}\bar{\rho}_{k}^{2}\left(x_{k}(h(\delta))-1\right)^{2},$

and by virtue of property B) and the proven convergence $h(\delta)\rightarrow
0$ the equality $\underset{\delta\rightarrow
0}{\lim}x_{k}(h(\delta))=x_{k}(0)=1$ is true. Therefore
$\tilde{A}_{h(\delta)}\rightarrow\bar{A}$. Now applying corollary 2 taking
into account (5), we obtain $\tilde{A}_{h(\delta)}^{+}\rightarrow\bar{A}^{+}$
and
$z_{\delta}=\tilde{A}_{h(\delta)}^{+}u_{\delta}\rightarrow\bar{z}=\bar{A}^{+}\bar{u}$
.$\square$

Now we study the question of improving the condition number of the original
matrix $\bar{A}$ when passing to the matrix $\tilde{A}_{h(\delta)}$.

###### Theorem 3

Let $r(\delta)$ be the rank of the matrix $\tilde{A}_{h(\delta)}$. Suppose
that for each $k$, $1\leq k\leq\bar{r}$, the relation $x_{k}(h)\sim 1+a_{k}h$
holds when $h\rightarrow+0$ , with $a_{1}<a_{2}<...<a_{\bar{r}}$. Then for
$\delta\rightarrow 0$ the following estimate is valid:
$\nu_{s}(\tilde{A}_{h(\delta)})\sim\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}(1-h(\delta)(a_{r(\delta)}-a_{1}))<\nu_{s}(\bar{A})$.

Proof. From the form of singular values of the matrix $\tilde{A}_{h(\delta)}$,
the equality
$\nu_{s}(\tilde{A}_{h(\delta})=\frac{\bar{)}{\rho}_{1}x_{1}(h(\delta))}{\bar{\rho}_{r(\delta)}x_{r(\delta)}(h(\delta))}$
follows, and $r(\delta)\leq\bar{r}$. Hence, taking into account the conditions
of the theorem and the convergence $h(\delta)\rightarrow 0$, we obtain a chain
of relations

$\nu_{s}(\tilde{A}_{h(\delta)})=\frac{\bar{\rho}_{1}x_{1}(h(\delta))}{\bar{\rho}_{r(\delta)}x_{r(\delta)}(h(\delta))}\sim\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}\frac{1+a_{1}h(\delta)}{1+a_{r(\delta)}h(\delta)}\sim\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}(1+a_{1}h(\delta))(1-a_{r(\delta)}h(\delta))\sim\\\
\sim\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}(1-h(\delta)(a_{r(\delta)}-a_{1}))<\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}\leq\frac{\bar{\rho}_{1}}{\bar{\rho}_{\bar{r}}}=\nu_{s}(\bar{A}),$

that proves the theorem.

Thus, the condition number of the matrix of the proposed method is better than
the condition number of the original exact matrix, at least for sufficiently
small $\delta$. If, in addition to the conditions of the theorem, the exact
matrix has a significant number of singular values close to zero, as in the
example 1, then it may turn out that
$\bar{\rho}_{r(\delta)}>>\bar{\rho}_{\bar{r}}$, and then
$\nu_{s}(\tilde{A}_{h(\delta)})<<\nu_{s}(\bar{A})$. This will be illustrated
below.

##  5 A specific method that implements the proposed scheme for solving SLAE

For simplicity, we assume that the singular values of the matrix $\bar{A}$ are
different: $\bar{\rho}_{k}>\bar{\rho}_{k+1},\leavevmode\nobreak\
k=1,2,...,M-1$.

As a central example, consider method using minimal pseudoinverse matrices
with condition number improvement or briefly, MPMI method. This method uses
the quantities

$x_{k}(h)=\left\\{\bar{x}_{k}(h),\leavevmode\nobreak\ 0\leq h\leq
h_{k};\leavevmode\nobreak\ 0,\leavevmode\nobreak\ h>h_{k}\right\\}$

for $1\leq k\leq\bar{r}$, where $\bar{x}_{k}(h)$ is the solution to the
equation
$\bar{x}_{k}^{4}-\bar{x}_{k}^{3}=h\bar{\rho}_{k}^{-4},\leavevmode\nobreak\
\bar{x}_{k}\in[1,\frac{3}{2}]$, $h_{k}=\frac{27}{16}\bar{\rho}_{k}^{4}$.

$\theta\left[x_{k}(h)\right]=\left\\{\frac{1}{\bar{x}_{k}(h)},\leavevmode\nobreak\
0\leq h\leq h_{k};\leavevmode\nobreak\ 0,\leavevmode\nobreak\
h>h_{k}\right\\}.$

Let us also define $H$: $H>h_{1}=\frac{27}{16}\bar{\rho}_{1}^{4}$. From the
definition of the function $x_{k}(h)$ and from the form of the function
$\bar{x}_{k}(h)$ the following properties are obtained: A) $1\leq
x_{k}(h)\leq\frac{3}{2}$, $1\leq h\leq H$; B)
$x_{k}(0)=x_{k}(+0)=\bar{x}_{k}(0)=1$; C) $x_{k}(h)$ is left continuous for
$0<h\leq H$. Property D) is also true. It follows from the increase in the
function $\bar{x}_{k}(h)$ as $h\in[0,h_{k}]$ (see [11, 12]). Then
$\theta\left[x_{k}(h)\right]=\frac{1}{\bar{x}_{k}(h)}>0$ decreases as
$h\in[0,h_{k}]$ and is equal to zero for $h>h_{k}$. Analyzing the asymptotic
behavior of the functions $x_{k}(h)$ for $h\rightarrow+0$, we can verify that
$x_{k}(h)\sim 1+h\bar{\rho}_{k}^{-4}$. This means that the conditions of the
theorem 3 are satisfied with $a_{k}=\bar{\rho}_{k}^{-4}$, and
$a_{k}=\bar{\rho}_{k}^{-4}<\bar{\rho}_{k+1}^{-4}=a_{k+1}$ for all admissible
$k$.

Thus, theorems 2 and 3 are valid, guaranteeing the convergence of the MPMI
method and the improvement of the condition number in it. In some cases, this
number differs significantly from $\nu_{s}(\bar{A})$. For example, the
theoretically possible case is
$x_{r(\delta)}(h(\delta))=x_{r(\delta)}(h_{r(\delta)})=\frac{3}{2}$ . It is
realized when the generalized solution of the equation
$\beta_{\delta}^{2}(h)=\delta^{2}+\mu_{\delta}^{2}$ is the discontinuity point
$h_{r(\delta)}$ functions $\beta_{\delta}(h)$. Then, as in theorem 3, we
obtain the relations

$\nu_{s}(\tilde{A}_{h(\delta)})=\frac{2}{3}\frac{\bar{\rho}_{1}x_{1}(h(\delta))}{\bar{\rho}_{r(\delta)}}\sim\frac{2}{3}\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}(1-h(\delta)(\bar{\rho}_{r(\delta)}^{-4}-\bar{\rho}_{k}^{-4}))<\frac{2}{3}\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}\leq\frac{2}{3}\nu_{s}(\bar{A}),$
(8)

and in this case, the condition number of the matrix used to solve the SLAE is
improved by at least one and a half times.

Next, we compare the MPMI method with some others.

TSVD method. Here, instead of $\bar{A}$, the matrix
$\tilde{A}_{r(\delta)}=U\tilde{R}_{r(\delta)}V^{T}$ is used, in which
$\tilde{R}_{r(\delta)}=\mathrm{diag}\left[\bar{\rho}_{1},\bar{\rho}_{2},...,\bar{\rho}_{r(\delta)},0,...,0\right]\leavevmode\nobreak\
$. Its rank $r(\delta)$ is found as a solution to the equation
$\beta^{2}_{\delta}(r)=\mu_{\delta}^{2}+\delta^{2}$ with monotonically non-
decreasing function $\beta^{2}_{\delta}(r)=\sum\limits_{k=r+1}^{m}v_{k}^{2}$.
We can verify that $r(\delta)=r$ for sufficiently small $\delta$. Thus, the
inequality
$\nu_{s}(\tilde{A}_{r(\delta)})=\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}\leq\frac{\bar{\rho}_{1}}{\bar{\rho}_{r}}=\nu_{s}(\bar{A})$
becomes equality for small $\delta$, and the TSVD method does not improve the
condition number.

Variants of the regularization method. In the Tikhonov regularization (TR), an
approximate solution of the SLAE is sought in the form
$z_{\delta}=\left(\alpha(\delta)I+\bar{A}^{\ast}\bar{A}\right)^{-1}\bar{A}^{\ast}u_{\delta}$,
where the parameter $\alpha(\delta)>0$ is selected by one of the known methods
(see, for example, [7, 12, 15, 16]). In any case, $\alpha(\delta)\rightarrow
0$ for $\delta\rightarrow 0$. Using the SVD for the matrix $\bar{A}$ we find
that $z_{\delta}=VT_{\delta}^{-1}U^{T}u_{\delta}$, where

$T_{\delta}=\mathrm{diag}\left[\frac{\alpha(\delta)+\bar{\rho}_{1}^{2}}{\bar{\rho}_{1}},\frac{\alpha(\delta)+\bar{\rho}_{2}^{2}}{\bar{\rho}_{2}},...,\frac{\alpha(\delta)+\bar{\rho}_{r}^{2}}{\bar{\rho}_{r}},0,...,0\right].$

Then, as in the proof of theorem 3, for $\alpha(\delta)\rightarrow 0$ we
obtain

$\nu_{s}(T_{\delta})=\frac{\bar{\rho}_{r}(\alpha(\delta)+\bar{\rho}_{1}^{2})}{\bar{\rho}_{1}(\alpha(\delta)+\bar{\rho}_{r}^{2})}=\frac{\bar{\rho}_{1}}{\bar{\rho}_{r}}\frac{(1+\alpha(\delta)/\bar{\rho}_{1}^{2})}{(1+\alpha(\delta)/\bar{\rho}_{r}^{2})}\sim\nu_{s}(\bar{A})\left(1+\frac{\alpha(\delta)}{\bar{\rho}_{1}^{2}}\right)\left(1-\frac{\alpha(\delta)}{\bar{\rho}_{r}^{2}}\right)\sim\\\
\sim\nu_{s}(\bar{A})\left(1-\alpha(\delta)\left(\frac{1}{\bar{\rho}_{r}^{2}}-\frac{1}{\bar{\rho}_{1}^{2}}\right)\right)<\nu_{s}(\bar{A}),$

taking into account that $\rho_{1}>\rho_{r}$.

Similar calculations cab be carried out for the method of solving SLAEs from
the work [16]:

$z_{\delta}=\left(\alpha(\delta)I+\bar{A}^{T}\bar{A}\right)^{-1}\bar{A}^{T}\bar{A}\bar{A}^{T}\left(\alpha(\delta)I+\bar{A}\bar{A}^{T}\right)^{-1}u_{\delta}=VM_{\delta}^{-1}U^{T}.$

As a result, we get a similar estimate:

$\nu_{s}(M_{\delta})\sim\nu_{s}(\bar{A})\left(1-\alpha(\delta)\left(\frac{1}{\bar{\rho}_{r}^{2}}-\frac{1}{\bar{\rho}_{1}^{2}}\right)\right)^{2}<\nu_{s}(\bar{A}).$

Thus, these variants of the regularization method theoretically provide, in
comparison with the original matrix $\bar{A}$, an improvement in the condition
number for the matrices $T_{\delta}$ and $M_{\delta}$, with the help of which
SLAEs are solved. However, these matrices have the same rank as
$\mathrm{rank}\bar{A}$, while in the MPMI method the rank $r(\delta)$ of the
matrix $\tilde{A}_{h(\delta)}$ can be significantly less than $\bar{r}$. This
possibility is demonstrated below in a numerical experiment. As a consequence,
the main terms of the condition numbers $\nu_{s}(T_{\delta})$ and
$\nu_{s}(M_{\delta})$, expressed by the term
$\frac{\bar{\rho}_{1}}{\bar{\rho}_{\bar{r}}}$ cannot be better than the
leading term of the condition number
$\nu_{s}(\tilde{A}_{h(\delta)})\sim\frac{\bar{\rho}_{1}}{\bar{\rho}_{r(\delta)}}$.
In addition, the estimate (8) for the number $\nu_{s}(\tilde{A}_{h(\delta)})$
may be fair. This guarantees a reduction of the condition number by one and a
half times. This feature is also demonstrated below.

##  6 Numerical experiments

Let us illustrate the main features of the proposed scheme for solving SLAE
using the example of the MPMI method discussed above. We solve a model SLAE
with an ill-conditioned matrix $\bar{A}$ of the form (4) and an exact solution
$\bar{z}=(1-x^{2})\sin(4\pi x),\,-1\leq x\leq 1$. The exact right-hand side of
this SLAE is calculated as $\bar{u}=\bar{A}\bar{z}$ and is perturbed by a
normally distributed random noise with zero mean so that the approximate
right-hand side $u_{\delta}$ satisfies the inequality
$\left\|u_{\delta}-\bar{u}\right\|_{2}\leq\delta\left\|\bar{u}\right\|_{2}$.
The number $\delta$ specified in this way allows one to estimate the error
level as a percentage. To solve the considered SLAE with different error
levels $\delta$, the MPMI method is used. For comparison, the TSVD and TR
methods were also used with the choice of the regularization parameter
$\alpha(\delta)$ according to the discrepancy principle [16].

Table 1 shows the results of calculations, namely the accuracy of the obtained
approximate solutions $z_{\delta}$:
$\Delta=\frac{\left\|z_{\delta}-\bar{z}\right\|_{2}}{\left\|\bar{z}\right\|_{2}}$
and condition numbers $\nu_{MPMI}=\nu_{s}(\tilde{A}_{h(\delta)})$,
$\nu_{TSVD}=\nu_{s}(\tilde{A}_{r(\delta)})$ and
$\nu_{TR}=\nu_{s}(T_{\delta})$.

Table 1. Comparison of solution accuracies and condition numbers

matrices of the MPMI, TSVD and TR methods  $\delta$ 0.005 0.01 0.05 0.1 0.2
0.3 $\Delta_{MPMI}$ 0.0024 0.0043 0.0117 0.0154 0.0333 0.0406 $\Delta_{TSVD}$
0.0027 0.0052 0.0131 0.0184 0.0346 0.0496 $\Delta_{TR}$ 0.0082 0.0108 0.0269
0.0358 0.0495 0.0989 $\nu_{MPMI}$ 20.972 20.971 10.353 10.353 10.353 5.6134
$\nu_{TSVD}$ 33.421 33.420 15.530 15.530 15.530 8.4172 $\nu_{TR}$ 2.3$\cdot
10^{12}$ 5.9$\cdot 10^{12}$ 3.3$\cdot 10^{13}$ 7.7$\cdot 10^{13}$ 2.6$\cdot
10^{14}$ 5.6$\cdot 10^{14}$

The table shows that the MPMI method has the best accuracy for all considered
levels $\delta$ of data disturbance. We especially note the high accuracy of
this method for large disturbances: $\delta=0.1,0.2,0.3$. Also, the table
clearly shows how much the condition numbers of the matrices of the MPMI
method are smaller than the condition numbers of other methods and,
especially, the condition numbers of the original SLAE matrix, $\nu_{s}\approx
3.37\cdot 10^{19}$. Note also that for the method from the work [16] the
results are very close to their analogues from the TR method and therefore are
not included in the table.

Calculations also showed that a significant decrease in the condition number
in the MPMI method (see (8)) occurs quite often. For example, for
$\delta=0.05$ it turns out that $h(\delta)$ is a discontinuity point of the
function $\beta_{\delta}(h)$, and therefore for such $h(\delta)$ the estimate
(8) is true (see fig.1B with
$\bar{\beta}_{\delta}(h)=\beta_{\delta}(h)/\|\bar{u}\|_{2}$). Effects of this
kind also occurred for other $\delta$.

In conclusion, we note that similar results were also obtained in numerical
experiments with other SLAEs that have ill-conditioned matrices.

This work was supported by the Russian Science Foundation (grant no.
23-41-00002).

## References

  * [1] Fox L. An introduction to numerical linear algebra, Oxford Univ. Press, NY, 1965.
  * [2] Forsythe G.E. and Moler C. Computer Solution of Linear Algebraic Systems, Prentice Hall, Englewood Cliffs, NJ, 1967.
  * [3] Lawson C.L. and Hanson R.J. Solving Least Squares Problems, Prentice Hall, Englewood Cliffs, NJ, 1974.
  * [4] Golub G.H. and Van Loan C.F. Matrix Computations, The John Hopkins University Press, 1989
  * [5] Voevodin Val.V. Numerical methods, Nauka, Moscow, 1977 (in Russian) .
  * [6] Tyrtyshnikov E.E. Methods of numerical analysis, Publishing center Academy, Moscow, 2007 (in Russian).
  * [7] Tikhonov A.N. and Arsenin V.Y. Solution of Ill-Posed Problems, Wiley, New York, 1977\.
  * [8] Varah J.M. On the numerical solution of ill-conditioned linear systems with applications to ill-posed problems, SIAM J. Num. Anal., V.10, 257267 (1973).
  * [9] Engl H.W., Hanke M. and Neubauer A. Regularization of Inverse Problems, Kluwer, Dordrecht, 1996.
  * [10] Golub G.H. Least squares, singular values and matrix approximation, Aplikace Matematiky, V.13, 4451 (1968).
  * [11] Leonov A.S. The minimal pseudo-inverse matrix method, U.S.S.R. Comput. Math. Math. Phys., V.27, No.4, 107117 (1987).
  * [12] Tikhonov A.N., Leonov A.S. and Yagola A.G. Nonlinear Ill-Posed Problems, V.2, Chapman and Hall, London, 1998.
  * [13] Leonov A.S. Solving Ill-Posed Inverse Problems: Outline of Theory, Practical Algorithms and Demonstrations in MATLAB, Librokom, Moscow, 2010 (in Russian).
  * [14] Albert A. Regression and the Moor-Penrose pseudoinverse, Academic Press, New York and London, 1972.
  * [15] Tikhonov A.N., Goncharsky A.V., Stepanov V.V. and Yagola A.G. Numerical Methods for the Solution of Ill-Posed Problems, Mathematics and Its Applications, Vol. 328, Springer, Dordrecht, 1995.
  * [16] Morozov V.A. Methods for Solving Incorrectly Posed Problems, Springer-Verlag, New York, 1984.


# IMDB-WIKI-SbS: An Evaluation Dataset for Crowdsourced Pairwise Comparisons

Nikita Pavlichenko
Yandex
Moscow, Russia
<EMAIL_ADDRESS>
Dmitry Ustalov
Yandex
Saint Petersburg, Russia
<EMAIL_ADDRESS>

###### Abstract

Today, comprehensive evaluation of large-scale machine learning models is
possible thanks to the open datasets produced using crowdsourcing, such as
SQuAD, MS COCO, ImageNet, SuperGLUE, etc. These datasets capture objective
responses, assuming the single correct answer, which does not allow to capture
the subjective human perception. In turn, pairwise comparison tasks, in which
one has to choose between only two options, allow taking peoples’ preferences
into account for very challenging artificial intelligence tasks, such as
information retrieval and recommender system evaluation. Unfortunately, the
available datasets are either small or proprietary, slowing down progress in
gathering better feedback from human users. In this paper, we present IMDB-
WIKI-SbS, a new large-scale dataset for evaluating pairwise
comparisons.111https://github.com/Toloka/IMDB-WIKI-SbS It contains 9,150
images appearing in 250,249 pairs annotated on a crowdsourcing platform. Our
dataset has balanced distributions of age and gender using the well-known
IMDB-WIKI dataset as ground truth. We describe how our dataset is built and
then compare several baseline methods, indicating its suitability for model
evaluation.

##  1 Introduction

Crowdsourcing has become an essential choice for producing large-scale
training and evaluation datasets for a broad spectrum of machine learning
tasks, including computer vision, natural language processing, etc. However,
crowdsourcing requires a careful task design that is usually reduced to a
classification task with a single known objective response (Zheng:17, ). In
such applications as information retrieval and recommender systems training
and evaluation, there is no objective response as the obtained rankings are
relatively complex to evaluate properly. This challenge is usually addressed
by relaxing the requirement of the objective answer and trying to gather the
_subjective human opinions_ on how the objects should be ranked (Beutel:19, ).

As most humans are non-experts in ranking evaluation, preferences are gathered
using pairwise comparisons (side-by-side, or SbS comparisons). In this
scenario, the crowdsourcing tasks offered to the workers show a pair of
objects with a question of which object is better. However, to train or
evaluate the systems described above, it is necessary to obtain the ranking
from the resulting comparisons. This process is called the _aggregation_ and
is fraught with various difficulties such as dealing with noisy or unreliable
responses from crowd workers. Scientific community has developed specific
models for aggregation of pairwise comparisons into ranked lists, including
Bradley-Terry (Bradley:52, ), CrowdBT (Chen:13, ), factorBT (Bugakova:19, ),
and others (Zhang:16, ). The resulting lists are then used as ground truth
datasets for computing quality criteria, such as normalized discounted
cumulative gain or expected reciprocal rate. Unfortunately, most evaluation
datasets reported in studies are either small but open (Janssens:12, ;
Chen:13, ; Zhang:16, ) or large but proprietary (Bugakova:19, ), which
complicates systematic studies of these algorithms and benchmarking models,
forcing the researchers to use non-realistic synthetic models in their
experiments.

Table 1: Summaries for related datasets and the dataset presented in this paper. Dataset | # of Objects | # of Pairs
---|---|---
PeopleAge (Janssens:12, ; Zhang:16, ) | 50 | 1,225
PeopleNum (Khan:14, ) | 39 | 741
EventTime (Zhang:16, ) | 100 | 4,950
ImageClarity (Zhang:16, ) | 100 | 4,950
Graded Readability (Chen:13, ) | 491 | 13,856
IMDB-WIKI-SbS [This Paper] | 9,150 | 250,249

We attempt in this paper to bridge the gap between these two extremes,
presenting IMDB-WIKI-SbS, a new large-scale dataset for evaluation pairwise
comparisons, building on the success of a well-known benchmark for computer
vision systems IMDB-WIKI (Rothe:18,
).222https://data.vision.ee.ethz.ch/cvl/rrothe/imdb-wiki/ Our dataset uses the
age information offered by IMDB-WIKI as ground truth while providing a
balanced distribution of ages and genders of people in photos. Table 1 shows
the comparison of our dataset to other datasets used in research papers in
terms of size.

The rest of the paper is organized as follows. Section 2 explains how IMDB-
WIKI-SbS was derived and sampled from the original IMDB-WIKI dataset. Section
3 describes how the sampled pairs of photos were annotated on the Toloka
crowdsourcing platform. Section 4 presents the empirical evaluation of three
baselines on our new dataset. Section 5 concludes with final remarks.

##  2 Dataset

This section describes how we chose the photos for our dataset from a popular
IMDB-WIKI dataset (Rothe:18, ) and produced sample pairs for further
crowdsourced annotation.

###  2.1 Choice of Photos

We have derived our dataset from the IMDB-WIKI as follows. As in the original
dataset, the odds ratio between males and females is roughly 1.4:1, and 50% of
the data represent people aged 28–45, which does not allow for unbiased
comparison between age and gender groups provided in the dataset. We decided
to narrow our dataset span from 10 to 70 years. We excluded all the entries in
which the face bounding box and gender label are missing. Then, we sampled
uniformly 75 photos per age and gender. One person can appear only once in
every group, but the same person can appear in multiple groups. If a person
has multiple photos in the group, we have chosen the one with the highest face
detector score provided by the dataset authors. As the result, our final
dataset contained $(60+1)\times 2\times 75=9,\\!150$ photos. We have noticed
that the original IMDB-WIKI contains the wrong gender information for some
photos. However, we believe the amount of such mistakes does not affect the
evaluation, given the size of our dataset.

###  2.2 Sampling Pairs

To produce a dataset allowing algorithm benchmarking, we decided to provide a
_connected_ graph of pairwise comparisons. To do so, we generated an Erdős-
Rényi graph containing 250,249 edges using the NetworkX toolkit (Hagberg:08,
). When generating the graph, we explicitly shuffled the pairs to ensure that
the distribution of labels would be uniform during crowd annotation. Although
pairwise comparisons graphs might be disconnected in real-world scenarios, we
believe that it is simpler to remove the available edges from our dataset and
then evaluate how well the algorithm recovers and handles the missing ones. We
decided to focus on the connected version due to the evaluation of the popular
Bradley-Terry model that requires such a graph configuration (Bradley:52, ).

##  3 Annotation

This section describes how we set up the quality control and annotated our
dataset on the Toloka crowdsourcing platform, https://toloka.ai/.

###  3.1 Quality Control

We used the following quality control configuration. We limited the minimum
time for task completion to at least four seconds. If the worker skipped more
than ten task pages, we blocked them from accessing additional tasks. Our
tasks were available to the top-50% workers on Toloka who set English as the
spoken language. Additionally, to avoid malicious behavior, we used 200 golden
comparisons that contained images of people with more than 20 years gap in
age. We mixed these golden comparisons with the regular tasks. Failing to
solve them correctly led to blocking access to our tasks.

###  3.2 Data Annotation

This dataset contains images of people with reliable ground-truth age assigned
to every image. The task interface contained two images and two buttons, left
and right. We intentionally omitted the “no answer” to avoid ambiguous cases
not handled by all the algorithms and prevent abusing this option by malicious
workers. This setup allowed us to obtain 250,249 comparisons by 4,091 workers.
One pair was annotated by only one unique worker, but a worker could annotate
multiple tasks. Each task page contained three regular tasks and one control
task in random order. While annotating, we used the cropped face images
instead of the complete photos. On average, each task page was done in 23
seconds, and there were about 20 task pages per worker. The whole annotation
process took 6 hours. As a result, 125,302 pairs labeled “left” and 124,947
pairs labeled “right”.

##  4 Evaluation

This section describes how we perform basic evaluation of baseline methods on
our dataset using the commonly used normalized discounted cumulative gain
(NDCG@k) score (Jarvelin:02, ).

###  4.1 Baselines

We have evaluated the performance of four baselines: a specific algorithm for
aggregation of pairwise comparisons, an algorithm that utilizes only
comparisons graph information, a random baseline, and the reversed ground-
truth ranking. We provide a brief description of these methods below.

#### Bradley-Terry.

Bradley-Terry is a probabilistic model for predicting the outcome of pairwise
comparisons (Bradley:52, ). Under this model, each item $i$ is associated with
some positive real-valued score $p_{i}$. Given two items $i$ and $j$ with
scores $p_{i}$ and $p_{j}$ respectively, probability of $i$ to be ranked
higher than $j$ is $\frac{p_{i}}{p_{i}+p_{j}}$. The likelihood of received
comparisons according to this probabilistic model is maximized with respect to
the scores. We use the implementation of the Bradley-Terry model available in
the Crowd-Kit computational quality control toolkit for crowdsourcing
(Ustalov:21, ).

#### PageRank.

PageRank is a classical algorithm for computing node centralities in a
directed graph (Brin:98, ). We built our directed graph as follows. As nodes,
we used all the 9,150 images in our dataset. As arcs, we used all the 250,249
pairs in our dataset. The direction of an arc represents which image in the
pair was labeled as older during the annotation. The resulting ranking is
obtained by sorting the items by their PageRank scores in descending order.

#### Random.

As a random baseline, we used a ranking obtained by sorting the items by a
uniformly generated random score. We ran this baseline to show the performance
gain obtained by specialized methods.

#### Reversed GT.

As a sanity check, we used a reversed ground truth (GT) dataset to ensure the
methods offer reasonable ranking order.

Table 2: Comparison of pairwise aggregation methods on the whole IMDB-WIKI-SbS dataset. Method | NDCG@100 | NDCG@100 (Female Only) | NDCG@100 (Male Only)
---|---|---|---
Bradley-Terry | $0.859$ | $0.879$ | $0.864$
PageRank | $0.803$ | $0.856$ | $0.864$
Random | $0.490$ | $0.533$ | $0.475$
Reversed GT | $\approx 0$ | $\approx 0$ | $\approx 0$

###  4.2 Results

Our evaluation in Table 2 shows that even the simple pairwise comparison
aggregation model performs well on our dataset compared to the random baseline
while leaving space for the more principled models. Additionally, we evaluated
the same baselines on two equally sized parts of our datasets separated by
gender and found no specific difference in algorithm performance. To assess
the difficulty of each age group, we evaluated the algorithm performance among
different age groups in Table 3. Even though each group contained the same
number of photos, aggregation algorithms performed worse than the random
baseline. We believe there are two reasons for this. First, for workers,
between-group comparisons were more straightforward to interpret than the in-
group ones, leading to the different grades of difficulty to the algorithms
and explaining the dramatic difference to the evaluation on the whole dataset.
Second, as the comparison graph, in this case, is not connected, the Bradley-
Terry model failed to deliver a meaningful ranking.

Table 3: Comparison on pairwise aggregation methods on different age groups in IMDB-WIKI-SbS. NDCG@10 is used due to the relatively small number of objects per group. Method | [10; 20] | [20; 30] | [30; 40] | [40; 50] | [50; 60] | [60; 70]
---|---|---|---|---|---|---
Bradley-Terry | $0.358$ | $0.461$ | $0.622$ | $0.436$ | $0.681$ | $0.570$
PageRank | $0.352$ | $0.682$ | $0.528$ | $0.629$ | $0.564$ | $0.715$
Random | $0.507$ | $0.501$ | $0.547$ | $0.517$ | $0.549$ | $0.577$

##  5 Conclusion

In this work, we have collected the largest available open dataset for
crowdsourced pairwise comparisons. The described IMDB-WIKI-SbS dataset has
been published at https://github.com/Toloka/IMDB-WIKI-SbS, and its labels are
made available under a CC BY license. We believe this dataset will foster the
development of better human-computer systems by allowing careful evaluation of
these systems against human judgments and encourage the scientific community
to develop better aggregation methods for crowdsourced pairwise comparisons.

## References

  * [1] Alex Beutel, Jilin Chen, Tulsee Doshi, Hai Qian, Li Wei, Yi Wu, Lukasz Heldt, Zhe Zhao, Lichan Hong, Ed H. Chi, and Cristos Goodrow.  Fairness in Recommendation Ranking through Pairwise Comparisons.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD ’19, page 2212–2220, Anchorage, AK, USA, 2019. Association for Computing Machinery.
  * [2] Ralph Allan Bradley and Milton E. Terry.  Rank Analysis of Incomplete Block Designs: I. The Method of Paired Comparisons.  Biometrika, 39(3/4):324–345, 1952.
  * [3] Sergey Brin and Lawrence Page.  The anatomy of a large-scale hypertextual Web search engine.  Computer Networks and ISDN Systems, 30(1):107–117, 1998.  Proceedings of the Seventh International World Wide Web Conference.
  * [4] Nadezhda Bugakova, Valentina Fedorova, Gleb Gusev, and Alexey Drutsa.  Aggregation of pairwise comparisons with reduction of biases.  In 2019 ICML Workshop on Human in the Loop Learning (HILL 2019), 2019.  arXiv:1906.03711 [cs.LG].
  * [5] Xi Chen, Paul N. Bennett, Kevyn Collins-Thompson, and Eric Horvitz.  Pairwise Ranking Aggregation in a Crowdsourced Setting.  In Proceedings of the Sixth ACM International Conference on Web Search and Data Mining, WSDM ’13, pages 193–202, Rome, Italy, 2013. Association for Computing Machinery.
  * [6] Aric A. Hagberg, Daniel A. Schult, and Pieter J. Swart.  Exploring Network Structure, Dynamics, and Function using NetworkX.  In Proceedings of the 7th Python in Science Conference, SciPy 2008, pages 11–15, Pasadena, CA, USA, 2008.
  * [7] Edouard Janssens.  1 to 100 Years Project - 100 Men, 2012.  https://edouardjanssens.com/art/1-to-100-years/men/.
  * [8] Kalervo Järvelin and Jaana Kekäläinen.  Cumulated Gain-Based Evaluation of IR Techniques.  ACM Transactions on Information Systems, 20(4):422–446, 2002.
  * [9] Asif R. Khan and Hector Garcia-Molina.  Hybrid Strategies for Finding the Max with the Crowd: Technical Report.  Technical report, Stanford University, 2014.
  * [10] Rasmus Rothe, Radu Timofte, and Luc Van Gool.  Deep Expectation of Real and Apparent Age from a Single Image Without Facial Landmarks.  International Journal of Computer Vision, 126(2):144–157, 2018\.
  * [11] Dmitry Ustalov, Nikita Pavlichenko, Vladimir Losev, Iulian Giliazev, and Evgeny Tulin.  A General-Purpose Crowdsourcing Computational Quality Control Toolkit for Python.  In The Ninth AAAI Conference on Human Computation and Crowdsourcing: Works-in-Progress and Demonstration Track, 2021.  https://www.humancomputation.com/assets/wips_demos/HCOMP_2021_paper_85.pdf.
  * [12] Xiaohang Zhang, Guoliang Li, and Jianhua Feng.  Crowdsourced Top-k Algorithms: An Experimental Evaluation.  Proceedings of the VLDB Endowment, 9(8):612–623, 2016.
  * [13] Yudian Zheng, Guoliang Li, Yuanbing Li, Caihua Shan, and Reynold Cheng.  Truth Inference in Crowdsourcing: Is the Problem Solved?  Proceedings of the VLDB Endowment, 10(5):541–552, 2017.