Nonlinear Kaczmarz algorithms and their convergence

Abstract

This paper proposes a class of randomized Kaczmarz algorithms for obtaining isolated solutions of large-scale well-posed or overdetermined nonlinear systems of equations. This type of algorithm improves the classic Newton method. Each iteration only needs to calculate one row of the Jacobian instead of the entire matrix, which greatly reduces the amount of calculation and storage. Therefore, these algorithms are called matrix-free algorithms. According to the different probability selection patterns of choosing a row of the Jacobian matrix, the nonlinear Kaczmarz (NK) algorithm, the nonlinear randomized Kaczmarz (NRK) algorithm and the nonlinear uniformly randomized Kaczmarz (NURK) algorithm are proposed. In addition, the NURK algorithm is similar to the stochastic gradient descent (SGD) algorithm in nonlinear optimization problems. The only difference is the choice of step size. In the case of noise-free data, theoretical analysis and the results of numerical based on the classical tangential cone conditions show that the algorithms proposed in this paper are superior to the SGD algorithm in terms of iterations and calculation time.

Introduction

Consider solving the following nonlinear system of equations

$$f\left(\mathit{x}\right)=\mathbf{0},$$

where $f:{\mathbb{R}}^n\to {\mathbb{R}}^m$ is a nonlinear vector-valued function, and $\mathit{x}\in {\mathbb{R}}^n$ is an unknown vector. Eqs. (1) can also be written in the form of component quantities:

$$f_i\left(\mathit{x}\right)=0,i=1,2,\dots ,m,$$

where at least one $f_i:\mathit{D}\left(f_i\right)\subseteq {\mathbb{R}}^n\to \mathbb{R}$ ($i=1,2,\dots ,m$) is a nonlinear operator. If ${\mathit{x}}^{\ast }\in {\mathbb{R}}^n$ exists such that $f\left({\mathit{x}}^{\ast }\right)=0$ holds, then ${\mathit{x}}^{\ast }$ is called the solution of the nonlinear system of Eqs. (1).

Nonlinear problems widely exist in many important fields such as national defense, economy, engineering, and management. For example, nonlinear finite element problems, nonlinear fracture problems, circuit problems, nonlinear programming problems can all be attributed to nonlinear problems [1], [2], [3]. As long as a differential equation contains the nonlinear term of the unknown function and its derivative, whether it is a difference method or a finite element method, the discretization will yield a nonlinear systems of equations. Furthermore, with the advent of the big data era, large-scale nonlinear problems are widely present in various application problems. Therefore, it is a common problem to solve large-scale nonlinear systems of equations of the form (1).

Assuming that the above equations have a unique solution in a certain region $\mathit{D}$ (i.e., the equations are consistent), the Newton method is the most commonly used method to solve this problem, but when the system is very large or the Jacobian matrix is singular, the Newton method will consume a lot of resources. At present, the algorithm for solving large-scale overdetermined nonlinear systems of equations is not perfect, but the Kaczmarz algorithm for solving large-scale overdetermined linear equations has been well developed. The classical Kaczmarz iteration composed of cyclic orthogonal projections was designed in 1937 by the Polish mathematician Stefan Kaczmarz [4]. The Kaczmarz algorithm mainly solves the following large-scale overdetermined consistent linear equations:

$$Ax=b,$$

where the matrix $A\in R^{m\times n}(m\ge n)$, $b\in R^m$, and $x$ is an $n$ dimensional unknown vector. The Kaczmarz method can be formulated as

$$x_{k+1}=x_k+\frac{b_i-〈a_i,x_k〉}{{\Vert a_i\Vert }_2^2}{a}_i,k=0,1,\dots$$

where $a_i$ denotes the $i$th row of $A$, and $b_i$ denotes the $i$th component of b, $i=k\left(modm\right)+1$.

Kaczmarz proved the convergence of the Kaczmarz method for (3) with an invertible coefficient matrix (m $=$ n) [4], but the convergence rate of the general Kaczmarz method is not easy to obtain, see [5], [6]. When the method sweeps through the rows of $A$ in a random manner rather than in the given order, Strohmer and Vershynin [7] proved this randomized Kaczmarz (RK) method converges with expected exponential rate when the row index $i$ is randomly chosen values in $\{1,2,\dots ,m\}$ with probability proportional to ${\Vert a_i\Vert }_2^2$.

Using the idea of the RK method, a nonlinear randomized orthogonal projection method—the NRK method is proposed. The iterative formula of the NRK algorithm and its origin are shown in Section 2. But in order to avoid calculating the entire Jacobian matrix, the NRK method does not directly use the RK algorithm’s random strategy of selecting projection rows, but uses $p_i=\frac{{\left|f_i\left(\mathit{x}\right)\right|}^2}{{\Vert f\left(\mathit{x}\right)\Vert }^2}$ as the probability criterion to select the projection row at each step, and then prove the linear convergence of the expected value. In addition, we are surprised to find that the numerical experiment with the use of the NK algorithm, which selects a row of the Jacobian matrix in order, is better, and the NK algorithm does not need to calculate the probability for each iteration, thereby reducing the amount of calculation and computing time (CPU). Similar to the NK algorithm, the NURK algorithm is proposed. The probability that the NURK algorithm selects a projection row each time is uniformly distributed. The linear convergence in expectation of the NURK algorithm is also proved.

The nonlinear successive over-relaxation (NSOR) algorithm for solving nonlinear optimization problems has been proposed as early as the 1960s, and its specific iterative formula can be seen in [8]. It actually solves the nonlinear system of equations which is derived from the problem of minimizing a strictly convex functional. In 1984, Brewster et al. presented the convergence analysis of the NSOR algorithm [9], [10]. In order to accelerate its convergence rate, Brewster et al. proposed different selection strategies of overrelaxation parameter $\omega$ [11], [12]. In addition, Ecker [13] and Vrahatis [14] et al. applied this method to solve nonlinear problems arised in various application. The NSOR algorithm is quite different from the NRK algorithm proposed in this paper, which can be compared in the following aspects. (1) Iteration formulas: The former only uses the diagonal elements of the Hessian matrix (i.e. the Jacobi matrix in this paper), and only changes one component of the unknown vector in each iteration. But the latter uses one row of Jacobi matrix to change all components of the unknown vector. (2) Scope of application: In the former method, the Hessian (Jacobian) matrix of a nonlinear problem is assumed to be positive or semi-positive definite, and its diagonal elements are not zero. However, the latter only needs Jacobian to be full column rank matrix. Furthermore, the former can only be used to solve the square nonlinear system of equations, but the latter can solve not only square but also overdetermined nonlinear system of equations. (3) Selection strategies of row index: The former method usually chooses the row index according to the order of the equations, while the latter method introduces a randomized selection strategy. In general, the NRK algorithm can be used widely, and it can achieve high convergence rate due to the introduction of randomized strategies.

Before this paper, Kowar et al. [15] introduced the idea of the Kaczmarz algorithm on the basis of Landweber iteration [16] and constructed the Landweber–Kaczmarz (LK) algorithm in 2007. Both the Landweber iteration and the LK algorithm are for solving nonlinear ill-posed problems. But the Landweber iteration is not suitable solving large-scale overdetermined nonlinear problems, and the LK algorithm’s step strategy is quite different from ours. The stochastic gradient descent (SGD) algorithm recently proposed by Bangti Jin et al. [17] can be regarded as a randomized version of the LK algorithm. The SGD algorithm selects the Jacobian rows according to the uniform probability each time, which is the same as the NURK algorithm. The difference between the two algorithms is that the NURK algorithm uses a exact step size for each iteration. However, even without the Assumption 2.1(ii)-(iv) in [17], the NURK algorithm can deliver the property of linear convergence. Moreover, Needell et al. also compared the randomized Kaczmarz algorithm with SGD algorithm in [18], and improved the RK algorithm to construct a randomized Kaczmarz algorithm with partially biased sampling. But if this algorithm is applied to the problem of solving nonlinear systems of equations, according to the idea of the nonlinear Kaczmarz algorithm mentioned later, the entire Jacobian matrix is needed, which will undoubtedly increase the amount of calculation. In order to avoid this problem, the NRK algorithm changes the probability of selecting rows of the Jacobian matrix, thereby speeding up the convergence speed.

A brief overview of this paper is given below: in Section 2, we introduce the NRK algorithm and give the proof of the convergence of the NRK algorithm; in Section 3, the NK and NURK algorithms and the corresponding proof of convergence are given; in Section 4, we compare this kind of nonlinear Kaczmarz algorithm with the SGD algorithm. in Section 5, we present numerical experiments using the three nonlinear Kaczmarz algorithms to solve classical nonlinear system of equations. The final Section 6 gives relevant summaries and conclusions.

The following notations are used in this paper. On the Euclidean space ${\mathbb{R}}^n$, let $\Vert \cdot \Vert$ be the Euclidean norm. If $\mathbf{A}\in {\mathbb{R}}^{m\times n}$, $\{{\mathit{a}}_1^T,\dots ,{\mathit{a}}_n^T\}$ denotes the set of the rows of $\mathbf{A}$. The $i$th row and $j$th column of matrix $\mathbf{A}$ is denoted by the $A_{ij}$  ($i,j\in \{1,2,\dots ,n\}$). $\Vert \mathbf{A}(i,:)\Vert$ $(i=1,2,\dots ,m)$ represents the norm of the $i$th row of matrix $\mathbf{A}$. ${\Vert \mathbf{A}\Vert }_F=\sqrt{{\sum }_{i,j}{A}_{ij}^2}$ means the Frobenius norm. $A^{+}$ is the generalized inverse of matrix $\mathbf{A}$. ${\Vert {\mathbf{A}}^{+}\Vert }_2=\frac{1}{{\sigma }_{min}}$, ${\sigma }_{min}$ is the smallest non-zero singular value of the matrix $\mathbf{A}$, ${\kappa }_F\left(\mathbf{A}\right)={\Vert \mathbf{A}\Vert }_F{\Vert {\mathbf{A}}^{+}\Vert }_2$.