Implicit iterative algorithm for solving regularized total least squares problems

Full Text

Introduction. Total least squares (TLS) are widely used in solving systems of linear algebraic equations with inaccurate data on the right and left sides.

Total least squares are widely used in many applied fields [1]. Including for system identification [2-5], image restoration [6, 7], tomography [8, 9], speech processing [10, 11].

There are many algorithms for solving total least squares problems. The classical algorithm for solving the total least squares problem based on SVD (singular value decomposition) [12]. The solution of the total least squares problem based on augmented systems is considered in [13, 14]. To solve large-scale linear systems of equations or linear systems of equations with a sparse matrix, iterative algorithms for total least squares are used: the Newton method [15, 16], Rayleigh iterations [17], Lanchotz iterations [18].

Various regularization methods are used to solve very ill-conditioned total least squares problems. Today, there are two main approaches to the regularization of total least squares problems: based on the truncated SVD [19] and Tikhonov’s regularization [20], as well as their modifications [21-24].

One way to improve the accuracy of the solution is to use iterative methods of regularization [25]. In [26], an implicit iterative algorithm for ordinary least squares based on SVD was proposed.

The condition number for total least squares is always greater than the condition number for ordinary least squares. Tikhonov’s regularization for total least squares makes it possible to approximate the condition number to the condition number of ordinary least squares [20].

This article proposes an implicit iterative algorithm to solve total least squares problems. When using the proposed algorithm, the condition numbers at each iteration turn out to be less than the condition numbers of ordinary least squares. This makes it possible to find the total least squares solution for very ill-conditioned problems.

It is proposed to use a restriction on the length of the solution vector as a stopping criterion for the iterative algorithm. The simulation results showed the high solution accuracy of the proposed implicit iterative algorithm to solve regularized total least squares problems.

1. Problem Statement. Let the overdetermined system of equations be defined as

$Ax=f\mathrm{,}$ (1)

where $A\in {ℝ}^{m\times n}$, $f\in {ℝ}^m$, $m>n$.

We will assume that the matrix $A$ and vector $f$ contain errors

$A=\tilde{A}+\Xi \mathrm{,}f=\tilde{f}+\xi \mathrm{.}$                                        

It is required to find a solution for the overdetermined system (1) using perturbed data $A$ and $f$.

To find an approximate solution vector from the errors data, the total least squares can be applied [12]. The total least squares approach minimizes the squares of errors in the values of both dependent and independent variables

$\underset{x}{\min }\parallel \begin{array}{cc}\xi & \Xi \end{array}{\parallel }_F\mathrm{,}\text{\hspace{0.33em}s.t.\hspace{0.33em}}\left(\tilde{A}+\Xi \right)x=\tilde{f}+\xi \mathrm{,}$                                

where $\left(\begin{array}{cc}\xi ,& \Xi \end{array}\right)$ is the augmented matrix and $\parallel \cdot {\parallel }_F$ is the Frobenius norm.

The solving of total least squares is reduced to finding the minimum of the objective function:

$\underset{x\in {ℝ}^n}{\min }\frac{\parallel Ax-f{\parallel }^2}{1+\parallel x{\parallel }^2}\mathrm{,}$ (2)

where $\parallel \cdot \parallel \parallel \cdot {\parallel }_2$ is the Euclidean norm.

The article proposes an implicit iterative algorithm for the regularized solution of the system of equations (1) according to the data with errors using the total least squares.

2. Iimplicit iterative algorithm for solving regularized TLS problems. Using the SVD, an arbitrary matrix  can be represented as follows:

$A=U\Sigma V^{\top }\mathrm{,}$ (3)

where $U=\left(\begin{array}{ccc}{u}_1& \cdots & u_n\end{array}\right)\in {ℝ}^{m\times n}$ and $V\left(v_1 ... v_n\right)\in {ℝ}^{n\times n}$ are orthogonal matrices; $\Sigma =\text{diag}\left(\begin{array}{ccc}{\sigma }_1\left(A\right)& \cdots & {\sigma }_n\left(A\right)\end{array}\right)$; ${\sigma }_1\left(A\right)\ge \cdots \ge {\sigma }_n\left(A\right)$ are singular numbers of matrix $A$; $u_i$ and $v_i$ are respectively left and right singular vectors of matrix $A$.

Let the augmented matrix of the system of equations be defined as

$\overline{A}=\left(A\mathrm{,}f\right)$.

A solution to the total least squares problem exists and is unique if the following condition is satisfied [27]:

$\sigma ={\sigma }_{n+1}\left(\overline{A}\right)<{\sigma }_n\left(A\right).$ (4)

When condition (4) is satisfied, the solution to problem (2) can be obtained from a biased normal system of equations [27]:

$\left(A^{\top }A-{\sigma }^2{E}_n\right)x=A^{\top }f\mathrm{.}$ (5)

Let $\mu$ be a positive constant. Equation (5) is equivalent to the following equation:

$\mu A^{\top }Ax+x=\mu {\sigma }^{2}x+x+\mu A^{\top }f\mathrm{.}$ (6)

The implicit iterative algorithm for equation (6) has the following form:

$\left({\mu }^{-1}{E}_n+A^{\top }A\right)x_{k+1}=\left({\sigma }^2+{\mu }^{-1}\right)x_k+A^{\top }f\mathrm{.}$ (7)

We write (7) as

$x_{k+1}={\left({\mu }^{-1}{E}_n+A^{\top }A\right)}^{-1}\left(\left({\sigma }^2+{\mu }^{-1}\right)x_k+A^{\top }f\right),$

or

$x_{k+1}={\Phi }_{\mu }{x}_k+g_{\mu }\mathrm{,}$ (8)

where ${\Phi }_{\mu }=\left({\sigma }_{\min }^2+{\mu }^{-1}\right){\left({\mu }^{-1}{E}_n+A^{\top }A\right)}^{-1}$, $g_{\mu }={\left({\mu }^{-1}{E}_n+A^{\top }A\right)}^{-1}{A}^{\top }f$.

Using the SVD decomposition of the matrix $A$ (3), let us perform the following transformations:

${\Phi }_{\mu }=\left({\sigma }^2+{\mu }^{-1}\right){\left({\mu }^{-1}{E}_n+A^{\top }A\right)}^{-1}=$

$=\left({\sigma }^2+{\mu }^{-1}\right)V{\left(\Sigma +{\mu }^{-1}{E}_n\right)}^{-1}{V}^{\top }=\left({\sigma }^2+{\mu }^{-1}\right)\sum _{i=1}^n\frac{v_i{v}_i^{\top }}{{\sigma }_i^2+{\mu }^{-1}}$;

$g_{\mu }={\left({\mu }^{-1}{E}_n+A^{\top }A\right)}^{-1}{A}^{\top }f=$

$={\left[V\left(\Sigma +{\mu }^{-1}{E}_n\right)V^{\top }\right]}^{-1}V\Sigma U^{\top }f=V{\left(\Sigma +{\mu }^{-1}{E}_n\right)}^{-1}{V}^{\top }V\Sigma U^{\top }f=$

$=V{\left(\Sigma +{\mu }^{-1}{E}_n\right)}^{-1}\Sigma U^{\top }f\underset{i=1}{\overset{n}{=\sum }}{v}_i\frac{{\sigma }_i}{{\sigma }_i^2+{\mu }^{-1}}{u}_i^{\top }f\mathrm{.}$

Then the implicit scheme (8) can be written based on the singular value decomposition in the following form:

$x_{k+1}=\left({\sigma }^2+{\mu }^{-1}\right)\sum _{i=1}^n\frac{v_i^{\top }{v}_i}{{\sigma }_i^2+{\mu }^{-1}}{x}_k+\sum _{i=1}^n\frac{{\sigma }_i{u}_i^{\top }{v}_i}{{\sigma }_i^2+{\mu }^{-1}}f\mathrm{,}k=0,1\dots \mathrm{.}$ (9)

3. Convergence and conditionality of an implicit iterative algorithm. The spectral radius of the transition matrix ${\Phi }_{\mu }$ is

$\rho \left({\Phi }_{\mu }\right)=\left(\mu {\sigma }^2+1\right){\lambda }_{\max }\left[{\left(E_n+\mu A^{T}A\right)}^{-1}\right]=\frac{\mu {\sigma }^2+1}{{\lambda }_{\min }\left(E_n+\mu A^{\top }A\right)}\frac{\mu {\sigma }^2+1}{1+\mu {\sigma }_n^2\left(A\right)}\mathrm{,}$

where ${\lambda }_{\max }$, ${\lambda }_{\min }$ are the maximum and minimum eigenvalues of the matrices, respectively.

The convergence condition of the implicit method of simple iterations (7) can be written as follows:

$\rho \left({\Phi }_{\mu }\right)=\frac{\mu {\sigma }^2+1}{1+\mu {\sigma }_n^{2}A}<1.$ (10)

If condition (4) is satisfied and $\mu >0$, condition (10) is always satisfied. This means that the iterative algorithm (8) converges for all cases where the biased normal system of equations has a unique solution.

It can be shown that the larger the value $\mu$, the higher the rate of convergence of the algorithm.

Let us show that algorithms (8) and (9) have different values of the condition numbers. The simple iteration method can be written as follows:

$x_{k+1}\underset{x\in {ℝ}^n}{=\mathrm{argmin}}{||\left(\begin{array}{c}A\\ \sqrt{{\mu }^{-1}}{E}_n\end{array}\right)x-\left(\begin{array}{c}f\\ \sqrt{{\mu }^{-1}+{\sigma }^2}{x}_k\end{array}\right)||}_2^2\mathrm{.}$ (11)

Formula (11) can be represented in the following form

$x_{k+1}=A_{\mu }^{+}{f}_{\mu }^{\left(k\right)}\mathrm{,}$

where $A_{\mu =}\left(\begin{array}{c}A\\ \sqrt{{\mu }^{-1}}{E}_n\end{array}\right)$, $f_{\mu }^{\left(k\right)}=\left(\begin{array}{c}f\\ \sqrt{{\mu }^{-1}+{\sigma }^2}{x}_k\end{array}\right)$; $A_{\mu }^{+}$ is a pseudoinverse Moore–Penrose matrix.

Since $\text{rank}\left(A_{\mu }\right)=n$, then $A_{\mu }^{+}$ can be calculated by the formula:

$A_{\mu }^{+}={\left(A_{\mu }^{\top }{A}_{\mu }\right)}^{-1}{A}_{\mu }^{\top }\mathrm{.}$                                          

In this case, the problem corresponds to the classical form of the implicit method of simple iterations:

$x_{k+1}={\left(A_{\mu }^{\top }{A}_{\mu }\right)}^{-1}{A}_{\mu }^{\top }{f}_{\mu }^{\left(k\right)}={\left(A^{\top }A+{\mu }^{-1}{E}_n\right)}^{-1}\left(A^{\top }f+\left({\mu }^{-1}+{\sigma }^2\right)x_k\right),$

${\kappa }_2\left(A^{\top }A+{\mu }^{-1}{E}_n\right)=\frac{{\lambda }_{\max }\left(A^{\top }A+{\mu }^{-1}{E}_n\right)}{{\lambda }_{\min }\left(A^{\top }A+{\mu }^{-1}{E}_n\right)}=\frac{{\sigma }_1^2+{\mu }^{-1}}{{\sigma }_n^2+{\mu }^{-1}}\mathrm{.}$                     

For the implicit method based on SVD decomposition, the condition number is equal to the condition number of the matrix $A_{\mu }$:

${\kappa }_2\left(A_{\mu }\right)={\left(\frac{{\sigma }_1^2+{\mu }^{-1}}{{\sigma }_n^2+{\mu }^{-1}}\right)}^{1/2}\mathrm{.}$

4. Stopping rule for an implicit iterative algorithm. There are a large number of stopping rules for iterative regularized algorithms [28-30]. In this article, we will use to stop the algorithm (5) the restriction on the value of the norm of the solution:

$\parallel x_{k+1}\parallel \le \delta \mathrm{,}$ (12)

where $\delta$ is the maximum allowable value of the Euclidean norm of the solution vector.

In contrast to Tikhonov’s total least squares regularization [20], condition (12) is verified directly without calculating indirect parameters.

5. Simulation results. Regularization Toolbox [31] was used to generate test cases. A matrix $A_{2000\times 4}$ with singular values $\sigma \left(\begin{array}{cccc}5\cdot {10}^{-4}& {10}^4& {10}^6& {10}^7\end{array}\right)$ was generated.

The true vector is $x_{\text{true}}={\left(\begin{array}{cccc}1& 1& 1& 1\end{array}\right)}^{\top }$.

The vector $f$ is $f=A_{2000\times 4}u$.

Gaussian white noise with zero mean and standard deviation ${\sigma }_f={\sigma }_A={10}^{-2}$ was added to the matrix $A_{2000\times 4}$ and the vector $f$.

The algorithm (5) was compared with the classical SVD-based TLS algorithm [12], the solution based on augmented systems [13], and regularized total least squares [20]:

$\left(A^{\top }A-{\sigma }^2{E}_n+\alpha E_n\right)x=A^{\top }f\mathrm{.}$ (13)

The condition number of the matrix $A^{\top }A-{\sigma }^2{E}_n+\alpha E_n$ is

${\kappa }_2\left(A^{\top }A-{\sigma }^2{E}_n+\alpha E_n\right)=\frac{{\sigma }_1^2-\left({\sigma }^2-\alpha \right)}{{\sigma }_n^2-\left({\sigma }^2-\alpha \right)}\mathrm{.}$

The parameter $\alpha$ was selected from the interval $\left(0,{\sigma }^2\right)$ with a step ${10}^{-4}{\sigma }^2$:

${\alpha }_i={10}^{-4}{\sigma }^{2}i\mathrm{,}i=0,1\dots \mathrm{,10000.}$                                 

The algorithms were compared by the relative mean square error (RMSE) of the solution

$\delta x_k=\frac{\parallel x_k-x_{\text{t}rue}{\parallel }_2}{\parallel x_{\text{t}rue}{\parallel }_2}\cdot 100\text{\hspace{0.05em}\%.}$

The simulation results are presented in Table 1. Figure 1 shows the relative root mean square error of the solution (8) in the $k$-th iteration for various values of the parameter ${\mu }^{-1}$. Figure 2 shows the relative root mean square error of solution (13) depending on the choice of parameter ${\alpha }_i$.

Table 1: RMSE of the solution

Figure 1: RMSE of the solution (8) at the k-th iteration for various values of the parameter ${\mathit{\mu }}^{\mathbf{-}\mathbf{1}}$: 1 – ${\mathit{\mu }}^{\mathbf{-}\mathbf{1}}\mathbf{=}{\mathbf{10}}^{\mathbf{-}\mathbf{5}}\mathit{\sigma }$, 2 – ${\mathit{\mu }}^{\mathbf{-}\mathbf{1}}\mathbf{=}{\mathbf{10}}^{\mathbf{-}\mathbf{5}}\mathit{\sigma }$; 3 –  ${\mathit{\mu }}^{\mathbf{-}\mathbf{1}}\mathbf{=}{\mathbf{10}}^{\mathbf{-}\mathbf{2}}\mathit{\sigma }$

Figure 2: RMSE of the solution (13) for various values of the parameter  ${\mathit{\alpha }}_{\mathbf{i}}\mathbf{=}{\mathbf{10}}^{\mathbf{-}\mathbf{4}}{\mathit{\sigma }}^{\mathbf{2}}\mathit{i}$

Conclusion. The paper proposes a new implicit iterative algorithm for solving regularized total least squares problems. The simulation showed that the proposed algorithm has a higher accuracy compared to the solutions obtained by total least squares algorithms, as well as the total least squares solution with Tikhonov regularization.

The proposed implicit iterative algorithm makes it possible to implement a constraint on the length of the solution vector without solving additional nonlinear equations.

The condition number of problems solved at each iteration is less than the condition number of systems with Tikhonov regularization.

Competing interests. We have no competing interests. Author’s ResponsibilitiesEach author has participated in the development of the concept of the article and in the writing of the manuscript. The authors are absolutely responsible for the submission of the final manuscript in print. Each author has approved the final version of the manuscript.

Funding. This work was supported by the Federal Agency of Railway Transport (projects nos. 122022200429-8, and 122022200432-8).

Acknowledgments. The authors thank the referees for careful reading of the paper and valuable suggestions and comments.

About the authors

Dmitriy V. Ivanov

Email: dvi85@list.ru
ORCID iD: 0000-0002-5021-5259
SPIN-code: 6672-4830
Scopus Author ID: 22937879800
http://www.mathnet.ru/person42123

Cand. Phys. & Math. Sci., Associate Professor, Dept. of Information Systems Security, Dept. of Mechatronics

Aleksandr I. Zhdanov

Author for correspondence.
Email: ZhdanovAleksan@yandex.ru
ORCID iD: 0000-0001-6082-9097
SPIN-code: 5056-3555
Scopus Author ID: 7102747969
ResearcherId: E-1433-2014
http://www.mathnet.ru/person41724

Dr. Phys. & Math. Sci., Proffessor, Dept. of Applied Mathematics & Computer Science

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Figure 1. RMSE of the solution (8) at the \(k\)-th iteration for various values of the parameter \(\mu ^{-1}\): 1 – \(\mu ^{-1}=10^{-5}\sigma \), 2 – \(\mu ^{-1}=10^{-1}\sigma \); 3 – \(\mu ^{-1}=10^{-2}\sigma\)

Download (85KB)

Indexing metadata

3. Figure 2. RMSE of the solution (13) for various values of the parameter \(\alpha _i=10^{-4}\sigma ^{2} i\)

Download (160KB)

Indexing metadata