ON NEWTON'S METHOD FOR SOLVING A SYSTEM OF NONLINEAR MATRIX EQUATIONS

Abstract

In this paper, we are concerned with the minimal positive solution to system of the nonlinear matrix equations $A_1X^2+B_1Y +C_1=0$ and $A_2Y^2+B_2X+C_2=0$, where $A_i$ is a positive matrix or a nonnegative irreducible matrix, $C_i$ is a nonnegative matrix and $-B_i$ is a nonsingular M-matrix for i = 1, 2. We apply Newton's method to system and present a modified Newton's iteration which is validated to be efficient in the numerical experiments. We prove that the sequences generated by the modified Newton's iteration converge to the minimal positive solution to system of nonlinear matrix equations.

1. Introduction

In this paper, we consider the following system of nonlinear matrix equations that can be expressed in the form

\(\left\{\begin{array}{l} F_{1}(X, Y)=A_{1} X^{2}+B_{1} Y+C_{1}=0 \\ F_{2}(X, Y)=A_{2} Y^{2}+B_{2} X+C_{2}=0 \end{array}\right.\)       (1.1)

where \(X, Y \in \mathbb{R}^{n \times n}\) are unknown matrices, \(A_i\) is a positive matrix or a nonnegative irreducible matrix, \(C_i\) is a nonnegative matrix for \(i\) = \(1,2, - B_1\) and \(-B_2\) are nonsingular \(M\)-matrices.

Set \(A=\left(\begin{array}{cc} A_{1} & 0 \\ 0 & A_{2} \end{array}\right), B=\left(\begin{array}{cc} B_{1} & 0 \\ 0 & B_{2} \end{array}\right), C=\left(\begin{array}{cc} C_{1} & 0 \\ 0 & C_{2} \end{array}\right)\), and \(Z= \left(\begin{array}{cc} X & 0 \\ 0 & Y \end{array}\right)\), then system (1.1) can be equivalently reformulated as

\(F(Z)=A Z^{2}+B P^{T} Z P+C=0\)       (1.2)

where \(P = \left(\begin{array}{cc} 0 & I_{n} \\ I_{n} & 0 \end{array}\right)\) is a permutation matrix.

In the past several decades, solving matrix equations has been a hot topic in linear algebra field. Among various matrix equations, the quadratic matrix equation

\(A Z^{2}+B Z+C=0\)       (1.3)

has draw much attention. This kind of matrix equation has important applications in Quasi-birth-death process, random walk, tandem Jackson queue, see \([1, 7, 11, 12] \)and the references there in.

The quadratic matrix equation of the form (1.2) is a generalization of equation (1.3), which has been widely studied. In \([2, 3]\), Newton’s method was applied to equation (1.3). Supporting theory and implication details were also showed there. Higham and Kim [5] incorporated exact line search into Newton’s method which improves the global convergence properties of Newton’s method. Moreover, there are many results about the perturbation analysis on equation (1.3), see \([5, 9, 10]\) for details about the normwise, mixed and componentwise condition numbers. Similar results can be obtained for the equation of the form (1.2).

An efficient way to solve system (1.1) is to apply Newton’s method to equation (1.2). Note that the size of the matrices in equation (1.2) is \(2 n \times 2 n\), the computational cost of Newton’s method may be very expensive if \(n\) is large. We present a modified Newton’s iteration where the matrices dealt with are of size \( n \times n\). We prove that the sequences generated by the modified iteration method converge to the minimal positive solution of system (1.1). It is validated in the numerical experiments that the modified Newton’s method works much more efficiently than the Newton’s method that is directly applied to equation (1.2).

This paper is organized as follows. In Section 2, we present a modified Newton’s method and prove that the sequences generated by the iteration converge to the minimal positive solution of equation (1.1). In Section 3, some numerical experiments are given to show the efficiency of the modified Newton’s iteration.

We begin with the notations used throughout this paper. \(\mathbb{R}^{n \times n}\) stands for the set of \(n \times n\) matrices with elements in field \(\mathbb{R}\). For \(X=\left(x_{i j}\right) \in \mathbb{R}^{n \times n}\), we write \(X \geq 0(X>0)\) and say that \(X\) is nonnegative (positive) if \(x_{i j} \geq 0\left(x_{i j}>0\right)\)holds for all \(i, j\), and \(X \geq Y(X>Y)\) is used as a different notation for \(X-Y \geq 0(X-Y>0).\) For a matrix \(A=\left(a_{i j}\right) \in \mathbb{R}^{n \times n}\), vec\((A)\) is a vector defined by vec\((A)\) \(=\left(a_{1}^{T}, \ldots, a_{n}^{T}\right)^{T}\)with \(a_i\) as the \(i\)-th column of A. For matrices \(B=\left(b_{i j}\right) \in \mathbb{R}^{m \times n}\) and \(C \in \mathbb{R}^{p \times q}, B \otimes C=\left(b_{i j} C\right)\) ) is the Kronecker product of matrices \(B\) and \(C\), which is a \(mp\) \(\times\) \(nq\) block matrix. In is the identity matrix of size \(n \times n\).

2. A modified Newton’s iteration

In this section, we propose a modified Newton’s method for obtaining the solution of system (1.1). We show that the sequences generated by the modified method converge to the minimal positive solution of (1.1).

Definition 2.1. . ([4]) \(A\) \(matrix\) \(A \in \mathbb{R}^{n \times n}\) is an \(M -matrix\) if \(A=s I-B\) for some nonnegative matrix \(B\) and \(s \geq \rho(B)\) where \(p\) is the spectral radius; it is \(a\) singular \(M -matrix\) if \(s = \rho(B)\) and \(a\) nonsingular \(M -matrix\) if  \(s > \rho(B)\).

Theorem 2.2. ([6]) The following are equivalent:

(1) \(A\) is a nonsingular \(M -matrix\).

(2) \(A^{-1}\) is nonnegative.

(3) \(Av\) \(> 0\) for some vector \(V > 0\),

(4) All eigenvalues of \(A\) have positive real parts.

(5) \(Av \geq 0\) implies \(v \geq 0.\).

One of the numerical methods to solve system (1.1) is to apply Newton’s method to equation (1.2). According to the ideas developed in \([2,3,5]\), Newton’s iteration for solving equation (1.2) can be stated as

\(\left\{\begin{array}{l} A Z_{i} H_{i}+A H_{i} Z_{i}+B P^{T} H_{i} P=-F\left(Z_{i}\right), \\ Z_{i+1}=Z_{i}+H_{i} \end{array} \quad i=1,2, \ldots\right.\)        (2.1)

Note that the matrices in equation (2.1) is of size \(2n \times 2n\), which implies that the computation cost by iteration (2.1) is very expensive if \(n\) is very large. Hence, it is not so practical. To reduce the computation cost, we propose a modified Newton’s iteration where the matrices dealt with are of size \(n \times n\)

Consider the following iteration

\(\left\{\begin{array}{l} A_{1} X_{i} H_{i, 1}+A_{1} H_{i, 1} X_{i}+B_{1} H_{i, 2}=-F_{1}\left(X_{i}, Y_{i}\right) \\ A_{2} Y_{i} H_{i, 2}+A_{2} H_{i, 2} Y_{i}+B_{2} H_{i, 1}=-F_{2}\left(X_{i}, Y_{i}\right) \\ X_{i+1}=X_{i}+H_{i, 1} \\ Y_{i+1}=Y_{i}+H_{i, 2} \end{array}\right.\)        (2.2)

Given \(X_0 = Y_0 = 0\), we will prove that the sequences \(\left\{X_{i}\right\}\) and \(\left\{Y_{i}\right\}\) generated by (2.2) converge to the minimal positive solution of system (1.1). By compatibility with Kronecker products, the first two equations of (2.2) are equivalent to

\(\begin{array}{l} \left(I_{n} \otimes A_{1} X_{i}+X_{i}^{T} \otimes A_{1}\right) \operatorname{vec}\left(H_{i, 1}\right)+\left(I_{n} \otimes B_{1}\right) \operatorname{vec}\left(H_{i, 2}\right)=\operatorname{vec}\left(-F_{1}\left(X_{1}, Y_{1}\right)\right) \\ \left(I_{n} \otimes A_{2} Y_{i}+Y_{i}^{T} \otimes A_{2}\right) \operatorname{vec}\left(H_{i, 2}\right)+\left(I_{n} \otimes B_{2}\right) \operatorname{vec}\left(H_{i, 1}\right)=\operatorname{vec}\left(-F_{2}\left(X_{1}, Y_{1}\right)\right) \end{array}\)        (2.3)

For convenience of notation, we let

\(M_{i}=-\left[\begin{array}{ccc} I_{n} \otimes B_{1} & I_{n} \otimes A_{1} X_{i}+X_{i}^{T} \otimes A_{1} \\ I_{n} \otimes A_{2} Y_{i}+Y_{i}^{T} \otimes A_{2} & I_{n} \otimes B_{2} \end{array}\right]\).

Then, (2.3) can be rewritten as

\(M_{i}\left[\begin{array}{l} \operatorname{vec}\left(H_{i, 2}\right) \\ \operatorname{vec}\left(H_{i, 1}\right) \end{array}\right]=\left[\begin{array}{c} \operatorname{vec}\left(F_{1}\left(X_{i}, Y_{i}\right)\right) \\ \operatorname{vec}\left(F_{2}\left(X_{i}, Y_{i}\right)\right) \end{array}\right]\)        (2.4).

Suppose that \(M_i\) is nonsingular, it follows from (2.2) that

\(\left[\begin{array}{c} \operatorname{vec}\left(Y_{i+1}\right) \\ \operatorname{vec}\left(X_{i+1}\right) \end{array}\right]=\left[\begin{array}{c} \operatorname{vec}\left(Y_{i}\right) \\ \operatorname{vec}\left(X_{i}\right) \end{array}\right]+\left(M_{i}\right)^{-1}\left[\begin{array}{c} \operatorname{vec}\left(F_{1}\left(X_{i}, Y_{i}\right)\right) \\ \operatorname{vec}\left(F_{2}\left(X_{i}, Y_{i}\right)\right) \end{array}\right]\),

which leads to

\(M_{i}\left[\begin{array}{c} \operatorname{vec}\left(Y_{i+1}\right) \\ \operatorname{vec}\left(X_{i+1}\right) \end{array}\right]=\left[\begin{array}{c} \operatorname{vec}\left(-A_{1} X_{i}^{2}+C_{1}\right) \\ \operatorname{vec}\left(-A_{2} Y_{i}^{2}+C_{2}\right) \end{array}\right]\)        (2.5)

Theorem 2.3. Suppose \(A_i\) is a positive matrix or a nonnegative irreducible matrix, \(C_i\) is a nonnegative matrix for \(i=1,2,-B_{1}\)  and  \(-B_2\) are nonsingular M-matrices in (1.1). Suppose that there is a pair of positive matrices \((U,V)\) such that \(F_{1}(U, V) \leq 0\) and \(F_{2}(U, V) \leq 0\). Set \(X_0 = Y_0 = 0\), then the sequences   \(\left\{X_{i}\right\}\) and \(\left\{Y_{i}\right\}\) generated by iteration (2.2) converge to the minimal positive solution of system (1.1), that is, there is a pair of matrices \((S_X,S_Y)\) which is the minimal positive solution of (1.1) such that

\(\begin{array}{l} \lim _{i \rightarrow \infty} X_{i}=S_{X} \\ \lim _{i \rightarrow \infty} Y_{i}=S_{Y} \end{array}\)

Moreover,

\(M_{i}=-\left[\begin{array}{cc} I_{n} \otimes B_{1} & I_{n} \otimes A_{1} X_{i}+X_{i}^{T} \otimes A_{1} \\ I_{n} \otimes A_{2} Y_{i}+Y_{i}^{T} \otimes A_{2} & I_{n} \otimes B_{2} \end{array}\right]\)

is a nonsingular M-matrix for each \(X_i,Y_i.\)

Proof. We use mathematical induction. Let U and V be positive matrices such that

\(\left\{\begin{array}{l} F_{1}(U, V)=A_{1} U^{2}+B_{1} V+C_{1} \leq 0 \\ F_{2}(U, V)=A_{2} V^{2}+B_{2} U+C_{2} \leq 0 \end{array}\right.\)         (2.6)

Since \(X_{0}=Y_{0}=0,\left(-B_{1}\right)^{-1} \geq 0,\left(-B_{2}\right)^{-1} \geq 0\) and \(H_{1}=P(-B)^{-1} C P^{T} \geq 0\) Therefore, the following statements

FIGURE 1. Comparision of CPU time

FIGURE 2. Comparision of CPU time and relative residue