# Solving Separable Nonlinear Equations Using LU Factorization

Authors: Yun-Qiu Shen; Tjalling J. Ypma
Journal: ISRN Mathematical Analysis
 (2013)

---

## Abstract

Separable nonlinear equations have the formF(y,z)≡A(y)z+b(y)=0, where the matrix A(y)∈Rm×N and the vector b(y)∈Rm are continuously differentiable functions of y∈Rn and z∈RN. We assume that m≥N+n, and F′(y,z) has full rank. We present a numerical method to compute the solution (y∗,z∗) for fully determined systems (m=N+n) and compatible overdetermined systems (m>N+n). Our method reduces the original system to a smaller system f(y)=0 of m−N≥n equations in y alone. The iterative process to solve the smaller system
only requires the LU factorization of one m×m matrix per step, and the convergence is quadratic. Once y∗ has been obtained, z∗ is computed by direct solution of a linear system. Details of the numerical implementation
are provided and several examples are presented.


---

## Body

## 1. Introduction

Many applications [1, 2] lead to a system of separable nonlinear equations:
(1)F(y,z)≡A(y)z+b(y)=0,
where the matrixA(y)∈Rm×Nand the vectorb(y)∈Rmare continuously differentiable functions ofy∈Rn andz∈RNwithm≥N+n. Typicallynis very small, and for compatible systems (i.e., those with exact solutions)mis usually close toN. We assume that F′(y,z)is Lipschitz continuous and has full rankN+nin a neighborhood of a solution(y*,z*); thus, we assume thatA(y)has full rankNin this neighborhood.Standard projection methods, such as VARPRO [3, 4], transform the problem (1) to the minimization of a function inyalone:
(2)miny∥[I-A(y)A+(y)]b(y)∥,
where we use the Euclidean norm throughout this paper.In [5, 6], we proposed a different method, using left orthonormal null vectors ofA(y), to reduce (1) to an equation of the form
(3)f(y)=0
inyalone. We assumed that the system is fully determined (m=N+n), andA(y)has full rankN. That algorithm was extended to overdetermined systems (m>N+n) with full rankA(y)in [7]. One QR factorization ofA(y)is required in each iterative step of the methods used to solve these smaller systems fory. For details of these methods and their relationship to other methods, see [1, 7].In this paper, we use a special set of linearly independent left null vectors ofA(y)to construct a bordered matrixM(y)=[A(y)∣K]which inherits the smoothness ofA(y). We use this to construct a system ofm-Nequations of the formf(y)=0in the n unknowns y alone. The smaller system inherits the Lipschitz continuity and nonsingularity of the Jacobian matrix of the original system (1) so that quadratic convergence of the Newton or Gauss-Newton method for solvingf(y)=0is guaranteed. The QR factorization ofA(y)used in previous methods is here replaced by an LU factorization ofM(y), so the cost of each iterative step may be significantly reduced. The method works for both fully determined systems and compatible overdetermined systems. This paper extends the work using bordered matrices and LU factorization for underdetermined systems in [8] (m<N+n) to fully determined and compatible overdetermined systems.In the next section, we show how to reduce the original system to a new system in the variablesyonly. The corresponding numerical algorithm is presented in Section 3. Examples and computational results are provided in Section 4.


## 2. Analysis

Assume that(y*,z*)is an exact solution of (1), and thatF′(y,z)is Lipschitz continuous and has full rankNin a neighborhood of(y*,z*). For allyin this neighborhood, we write the singular value decomposition ofA(y)as
(4)A(y)=U(y)Σ(y)VT(y)≡[u1(y),…,um(y)][σ1(y)⋱σN(y)0⋯0⋮0⋯0]×[v1(y),…,vN(y)]T,
whereU(y)andV(y)are orthogonal matrices andσ1(y)≥⋯≥σN(y)>0. Since the singular valuesσ1(y),…,σN(y)are continuously dependent ony, andA(y*)has full rankNand thusσN(y*)>0, we have
(5)σN(y)≥12σN(y*)
in a small neighborhood ofy*. In this neighborhood, the following result holds.Theorem 1.
LetA(y) be continuously differentiable and full rank in a neighborhood of y*. Then, there exists an m×(m-N) constant matrixKsuch that the bordered matrix
(6)M(y)=[A(y)∣K]
has full rank in a neighborhood ofy*.Proof.
Lety-be a point neary*for which the singular valueσN(y-)satisfies (5). Let
(7)K=[uN+1(y-),…,um(y-)].
Then,(8)M(y-)=[A(y-)∣K]=[U(y-)Σ(y-)VT(y-)∣uN+1(y-),…,um(y-)]=U(y-)[[σ1(y-)⋱σN(y-)]VT(y-)0N×(m-N)0(m-N)×NI(m-N)×(m-N)]
which implies that M(y-)is invertible with∥M-1(y-)∥≤max{1,2/σN(y*)}by (5). For ally,y-in a sufficiently small neighborhood ofy*continuity ensures that we can satisfy
(9)∥M(y)-M(y-)∥=∥[A(y)-A(y-)∣0m×(m-N)]∥<ϵ=1(max{1,2/σN(y*)}).
This guarantees that for all suchy,y-(10)∥[M(y-)]-1(M(y-)-M(y))∥≤(max{1,2σN(y*)})ϵ<1;
hence by the Neumann Lemma [9],M(y)is invertible in this neighborhood.SinceM(y)in invertible, there exists anm×NmatrixP(y)and anm×(m-N)matrixC(y)with linearly independent columns that satisfy
(11)[PT(y)CT(y)]M(y)=Im.We use this to reduce the original system (1) to the form (3) for fully determined systems and compatible overdetermined systems with full rankF′(y,z) as follows.Theorem 2.
(a)F(y,z)=A(y)z+b(y)=0if and only if
(12)f(y)≡CT(y)b(y)=0,z=-ζ(y)≡-PT(y)b(y).
(b)F′(y*,z*)has full rank if and only iff′(y*)has full rank.Proof.
(a) Using (6) and (11), we have
(13)[PT(y)CT(y)][A(y)z+b(y)]=[z+PT(y)b(y)CT(y)b(y)],
which implies the result since[PT(y)CT(y)]is nonsingular.
(b) Differentiating (13) and usingF(y*,z*)=A(y*)z*+b(y*)=0, we have
(14)[PT(y*)CT(y*)]F′(y*,z*)=[ζ′(y*)INf′(y*)0(m-N)×N],
which implies thatF′(y*,z*)has full rank if and only iff′(y*)has full rank.


## 3. Computation

Based on (12), we solve them-Nequations innunknownsf(y)=0fory*and then computez*=-ζ(y*)=-PT(y*)b(y*). For the fully determined system, we use Newton’s method
(15)f′(y(k))(y(k+1)-y(k))=-f(y(k))
to computey*, and for compatible overdetermined systems, we use the Gauss-Newton method
(16)[(f′(y(k)))Tf′(y(k))](y(k+1)-y(k))=-(f′(y(k)))Tf(y(k)).To evaluate(17)f(y)=CT(y)b(y),
we use an LU factorization ofM(y)≡[A(y)∣K]to solve forCT(y)from (11):
(18)CT(y)M(y)=[0(m-N)×N∣Im-N].Now, considerf′(y). Use the following notation:(B(y))j′denotes a matrix of the same size as a given matrixB(y)but with each entry being the partial derivative of the corresponding entry ofB(y)with respect toyj. Then,
(19)f′(y)=[(CT(y))1′b(y),…,(CT(y))n′b(y)]+CT(y)b′(y).Differentiating (18), the matrices(CT(y))j′,j=1,…,ncan be computed by solving
(20)(CT(y))j′M(y)=-CT(y)[(A(y))j′∣0],
using the previously computedCT(y)and the LU factorization ofM(y).Notice that Lipschitz continuity ofF′(y,z)implies Lipschitz continuity of(A(y))j′,j=1,…,nandb′(y), and also ofM(y)andM-1(y). By (18) this implies Lipschitz continuity ofCT(y)and (20) then implies Lipschitz continuity of(CT(y))j′,j=1,…,n. Using (19), it follows thatf′(y)is Lipschitz continuous. Now, Theorem 2(b) guarantees quadratic convergence of the Newton and Gauss-Newton methods fory(0)sufficiently neary*since the corresponding convergence conditions are satisfied [9–12]. In particular, for the Gauss-Newton method, our assumption that we have a compatible overdetermined system implies that we have a zero residual problem.An outline of our algorithm follows.Algorithm 3.
Given:F(y,z)=A(y)z+b(y)with the corresponding positive integersm,N,n, a small positive numberϵ, and a pointy(0)near the solutiony*. Compute an SVD factorization ofA(y(0))and use it to formKas in (7). Fork=0,1,2,…, do steps (a)–(e).(a)
FormM(y(k))as in (6) and compute its LU factors.(b)
Formf(y(k))as in (17) and (18).(c)
Formf′(y(k))as in (19) and (20).(d)
Findy(k+1)using (15) ifm=N+nor (16) ifm>N+n.(e)
If∥y(k+1)-y(k)∥<ϵstop, outputy*≈y(k+1)and obtainz*from (12). Otherwise replacekbyk+1and go to(a).
By (11), in step(e)z*=-ζ(yk+1)=-PT(yk+1)b(yk+1)=-[IN,0]M-1(yk+1)b(yk+1). Thus, z*is the firstNcomponent of-M-1(yk+1)b(yk+1), which can be computed by using an LU decomposition ofM(yk+1).
Each iteration of our method requires one LU factorization ofM(y(k))or about(2/3)m3flops [13]. A QR factorization of them×NmatrixA(y(k))costs about2N2(m-N/3)flops [13], so ifmis close toNthe new method approximately halves this cost. The matrix to be factored in step(d)is onlyn×nand typicallynis very small, so this cost is negligible. SincemandNare typically large, the computational cost of the other steps is also small relative to the cost of the matrix factorization.


## 4. Examples

The following three examples illustrate the method.Example 1.
Consider(21)A(y)z+b(y)≡[y12+y22-y1-y2-2222y12+y2200y1+y2]z+[00-2y1+y21]=0.
This equation has a solution(y1*,y2*)=(2,-1), (z1*,z2*)=(1,-1), at whichA(y*)has full rankN=2. We choosey-=y(0)=[2.1-1.1]. An SVD ofA(y(0))produces the matrix block
(22)K=[0.2935-0.67800.16680.7264-0.19620.0576-0.9206-0.0968].We list the errors iny(k)in Table 1, which clearly shows quadratic convergence of the Newton iteration. Usingy(3)≈y*givesz*≈-ζ(y(3))with||-ζ(y(3))-z*||=4.0868e-12.Table 1
Newton iterations for Example1.
k
y
1
(
k
)
-
y
1
*
y
2
(
k
)
-
y
2
*
∥
y
(
k
)
-
y
(
k
-
1
)
∥
0
1.0000
e
-
001
-
1.0000
e
-
001
—
1
-
2.3657
e
-
003
4.4412
e
-
003
1.4624
e
-
001
2
-
5.7794
e
-
006
1.8337
e
-
008
5.0292
e
-
003
3
2.5508
e
-
012
7.5073
e
-
012
5.7794
e
-
006The next example is a discretized version of an interface eigenvalue problem for a differential equation on two intervals. For references on interface problems in differential equations, see [14].Example 2.
Consider(23)A(y)z+b(y)≡[y1E2k+1-I2k+10(2k+1)×(2k+1)0(2k+1)×(2k+1)y2E2k+1-I2k+1ek+1T00ek+1T]z+[0(2k+1)×10(2k+1)×1-11]=0,
where the(2k+1)×(2k+1)tridiagonal matrix
(24)E2k+1=[2-1-12-1…-12].
For eigenvalues of such matrices, see [15]. We usek=9so thatn=2,m=40, andN=38. This equation has a solution (y*)T=(csc2(π/40)/4,csc2(3π/40)/4) and fori=1,…,19,z*(i)=sin(iπ/20)andz*(i+19)=sin(3iπ/20), at whichA(y*)has full rankN=38. We choosey-=y(0)=[41.54.75]. An SVD ofA(y(0))produces an appropriate40×2matrix blockK. We list the errors iny(k)in Table 2, again showing quadratic convergence. Usingy(3)≈y*givesz*≈-ζ(y(3))with ∥-ζ(y(3))-z*∥=9.9027e-15.Table 2
Newton iterations for Example2.
k
y
1
(
k
)
-
y
1
*
y
2
(
k
)
-
y
2
*
∥
y
(
k
)
-
y
(
k
-
1
)
∥
0
8.8809
e
-
001
1.6257
e
-
001
—
1
-
4.2544
e
-
003
-
2.3859
e
-
003
9.0746
e
-
001
2
-
9.7644
e
-
008
-
5.1531
e
-
007
4.8774
e
-
003
3
7.1054
e
-
015
-
2.3981
e
-
014
5.2448
e
-
007The relative flop count(25)[2N2(m-N/3)-2m3/3][2N2(m-N/3)]
shows that in this example using one LU factorization ofM(y(k))instead of a QR factorization ofA(y(k)), as used in alternative methods, reduces the flop count per iteration by approximately46%. In general, ifmandNare large and close to each other, as is usually the case in compatible systems, the flop count is approximately halved.We now give an example of a compatible overdetermined system.Example 3.
Consider(26)A(y)z+b(y)≡[A1(y1,y2)A2(y1,y2)]z+[b1(y1,y2)b2(y1,y2)]=0,
where the(N-2)×NmatrixA1(y1,y2), the5×NmatrixA2(y1,y2), andb2(y1,y2)∈R5are, respectively,
(27)[y1-y21y1-y21⋱⋱⋱y1-y21],[0⋯0y121002y20⋯01100-y1y220011],[-5-2-y2-2-7-y1]
andb1(y1,y2)=0∈RN-2. Here, n=2andm=N+3.We chooseN=77so thatm=80. This equation has a solution(y1*,y2*)=(2,3)andzj*=1forj=1,…,77, at whichA(y*)has full rank. We choosey-=y(0)=[2.23.2]. An SVD ofA(y(0))produces an appropriate80×3matrix blockK. We list the errors iny(k)in Table 3, clearly showing quadratic convergence of the Gauss-Newton method. Usingy(4)≈y*givesz*≈-ζ(y(4))with∥-ζ(y(4))-z*∥=3.4101e-13. In each iteration of this example, using one LU factorization ofM(y(k))instead of one QR factorization ofA(y(k))reduces the flop count by approximately47%, measured by (25).Table 3
Gauss-Newton iterations for Example3.
k
y
1
(
k
)
-
y
1
*
y
2
(
k
)
-
y
2
*
∥
y
(
k
)
-
y
(
k
-
1
)
∥
0
2.0000
e
-
001
2.0000
e
-
001
—
1
-
6.7103
e
-
003
-
1.7935
e
-
004
2.8775
e
-
001
2
-
6.2985
e
-
006
-
1.7406
e
-
005
6.7060
e
-
003
3
-
3.5553
e
-
011
6.2474
e
-
011
1.8511
e
-
005
4
0.0000
e
+
000
-
8.8818
e
-
016
7.1883
e
-
011


## 5. Conclusions

We present an algorithm for solving the separable equation(28)F(y,z)≡A(y)z+b(y)=0,
where the matrixA(y)∈Rm×Nand the vectorb(y)∈Rmare continuously differentiable functions ofy∈Rn andz∈RNwithm≥N+n. We assume that F′(y,z)is Lipschitz continuous and has full rankN+nin a neighborhood of a solution(y*,z*); thus, A(y)has full rankNin this neighborhood. Our technique replaces (28) bym-Nequations innunknownsf(y)=0which we solve by either the Newton or the Gauss-Newton iterative method, in both cases with quadratic convergence. Our method uses one LU factorization of anm×mmatrix
(29)M(y)=[A(y)∣K]
instead of a QR factorization of the matrixA(y)per iteration, and may thus be substantially more efficient than competing methods (approximately halving the cost) whenmis close toNandNis large, as is usually the case in compatible systems. The method is applicable to fully determined separable systems and to compatible overdetermined separable systems.



---

*Source: 258072-2013-06-24.xml*