Solution of Nonlinear Least-Square Problems
on a Multiprocessor
Thomas F. Coleman1
Paul E Plassmann2
TR 88-923
June 1988
Department of Computer Science
Cornell university
Ithaca, NY 14853-7501
1 Research partially supported by the Applied Mathematical Sciences Research Program (KC-o4-o2)
of the Office of Energy Research of the u.s Department of Energy under grant DE-FGo2-86ER25013.A
000.
2 Research partially supported by the Gomputational Mathematics Program of the National Science
Foundation under grant DMS-8706133.
Solution of Nonlinear Least-Squares Problems on a Multiprocessor
Thomas F. Coleman1
Paul E. Plassmann2
Computer Science Department &
Center for Applied Mathematics
Cornell University
Ithaca, New York 14853
June 1988
Abstract: In this paper we describe algorithms for solving nonlinear least-squares
problems on a message-passing multiprocessor. We demonstrate new parallel algo-
rithms, including an efficient parallel algorithm for determining the Levenberg-
Marquardt parameter and a new row-oriented QR factorization algorithm. Experimen-
tal results obtained on an Intel iPSC hypercube are presented and compared with
sequential MINPACK code executed on a single processor. These experimental results
show that essentially full efficiency is obtained for problems where the row size is
sufficiently larger than the number of processors. These algorithms have the advantage
of involving only simple data movements and consequently are not constrained to the
hypercube architecture.
Key Words: hypercube
message-passing multiprocessor,
algorithms
computer, Levenberg-Marquardt, nonlinear least-squares,
parallel algorithms, QR-factorization, trust-region
Subject Classification: AMS/MOS: 65H10, 65F05, 65K05, 65K1O, 9OC3O
Research partially supported by the Applied Mathematical Sciences Research Program (KC-04-02) of
the office of Energy Research of the U.S. Department of Energy under grant DE-FGO2-86ER25013.AOOO.
2 Research partially supported by the Computational Mathematics program of the National Science
Foundation under grant DMS-8706l33.
1. Introduction
A common computational problem is the minimization of the function `v: R ? R where
`v is the sum of squares of nonlinear functions. That is, `v can be expressed in terms of a func-
tionF :R? ?R?, m>n, by the equation:
`v(x) = `/2 II F(x) 1122 = `/2 ?			(1.1)
i=1
where f? is the I -th component of F. In this paper we will describe a parallel implementation
of the Levenberg-Marquardt algorithin [M77] for solving these nonlinear least-squares prob-
lems. The experimental results presented in this paper were obtained on a hypercube multipro-
cessor but the algorithms themselves are more general. In fact, all that is required of multipro-
cessor interconnection topology is support of a ring embedding and means for efficient gather
and broadcast operations.
Unlike in the parallel solution of systems of nonlinear equations [CL87, BSS88], low
rank updates to the Jacobian approximation are usually not used in the small residual nonlinear
least-squares setting [GMW8l]. Consequently, the dominant costs in these problems are the
approximation of the Jacobian J (xk), where Xk is the k -th iterative approximation to the solu-
tion, and the computation of the QR factorization of this approximation to J(Xk).
The development of parallel algorithms for these two aspects of the problem, the approxi-
mation and QR factorization of J, represents a dilemma when considering how to distribute the
elements of the Jacobian onto the processors. This dilemma arises because the approximation
of J by finite differences is most naturally approached by a column-oriented method and the
QR factorization stage suggests a row-oriented solution. In a column-oriented algorithm for
approximating the Jacobian the columns of J are partitioned into sets that assign computation
of columns to processors by these sets. Assuming that F (Xk) is globally available, the j -th
column of J can be approximated by computing F(x?+Tej), where ej is the j-th column of the
identity matrix and T a constant, for each column I assigned to a processor. These computa-
tions are entirely local and if n/p?l, where p is the number of processors, then the computa-
tional load can be well-balanced by an even distribution of work to the processors. On the
other hand, because we have that m >??, we expect a row-oriented approach, in which the rows
of J are assigned to processors, to outperform a column-oriented method in the QR factoriza-
tion stage. Using a row-oriented method in this case results in an algorithin whose efficiency
-2-
depends on the ratio mIp rather than n Ip.
We have chosen to pursue a row-oriented algoritlim because experience has shown that
computational costs involved in the QR factorization stage often dominate the Jacobian approx-
imation stage. In addition, it is often the case that evaluation of the function F is separable.
`mat is, if !? is the set of row indices assigned to processor i then the evaluation of F (x) can
be effectively broken up into blocks, F?j(X ) = (fi (x) j E j?), where evaluation of each block is
sufficienfly independent to allow for parallel evaluatio?
The distribution of the rows of the Jacobian onto the processors determines the basic
communication structure of the algorithins used in an implementation. We have chosen to use
a ring embedding (e.g. a Gray code ordering of the nodes on a hypercube) and wrap the rows
of the Jacobian onto this ring. Specifically, if the processors on the ring are numbered
O,l,2;...,p-l then row k of the Jacobian would be assigned to processor (k--Hl)mod(p). Note,
however, that components of F may be reordered to facilitate separability. This observation is
especially applicable in the sparse setting where grouping rows with similar norizero structure
on the same processor would help ensure separable evaluation. Such an assignment would also
make sense for a sparse QR factorization since much progress in the triangular reduction of the
system could be made by locally applied Givens rotations.
The remainder of this paper is organized as follows. In section 2 we briefly review the
relevant aspects of the Levenberg-Marquardt algorithm necessary to explain details of our algo-
rithins. In section 3 we consider the problem of approximating the Jacobian and in section 4
we present a new row-oriented parallel QR factorization. Section 5 describes our parallel algo-
rithin for determining the Levenberg-Marquardt parameter. Finally, we present experimental
results and conclusions in section 6.
-3-
2. Basics of the Levenberg-Marquardt Algorithm
In this section we consider the essential aspects of the Levenberg-Marquardt algorithm
relevant to a parallel perspective. For a more detailed description of the algorittirn we refer the
interested reader to the excellent article by Mor? [M78]. As stated previously, the nonlinear
least-squares problem is to minimize ?(x) as given in equation (1.1). Assuming that each
f? E C2, then the gradient and Hessian of'j, are given by the expressions:
= JTF			(2.1)
and
v2'II=JTJ+?hv2L. ,			(2.2)
i=1
where J is the Jacobian of F. When one expects the value of the ftmction to he small at the
solution, known as small residual problems [MGW8 1 j, the Hessian can be approximated by
jTj in a neighborbood of the solution. This approximation, which ignores the second term in
the expression for the Hessian, represents a significant computational savings.
Based on a quadratic model to the function `ir using the above Hessian approximation the
Levenberg-Marquardt algorithin [M77, C84] solves a sequence of tmst-region problems of the
form
min ( iiJp +F 112 St. lIP 112 <?A )			(2.3)
A solution to this problem is obtained in two steps. First, the Gauss-Newton step, PGN is
determined by solving the problem
min ( Jp + F 112 )			(2.4)
by computing a QR factorization of j and solving the upper triangular least-squares			system
[R0]p = QTF			(2.3)
If IIPGN II <--H A, then the step is accepted,1 otherwise we determine the Lagrange multiplier ?
and vector p*, lIp* II = A, such that
(JTJ+?I)Th =?JTF
1 contingen? of course, upon sufficient decrease in ?
(2.6)
-4-
Using the previously computed QR factorization of J and given a value of ?, the solution of
the least-squares system
R
0 p(? = Q0TF
xiAi
(2.7)
yields the value of the function ?? = II? (? I --H A. Since ? is monotonically decreasing on
[O,oo) from the positive value of IIPGN --H A at ?=O to the negative number --HA as ??oo, an
appropriate implementation of Newtons method is guaranteed to find the zero of ? at ?, the
Lagrange multiplier [M78, MGH8O].
From this discussion we note that there are three main computational tasks that need to
be addressed in a parallel implementation of the Levenberg-Marquardt algorithm. These
features are: (1) the approximation of the Jacobian J (x), (2) the QR factorization of J (x)
necessary to solve equation (2.5), and (3) computation of the Levenberg-Marquardt parameter
that solves equation (2.6). In the following sections we present our implementation ot tne first
task and new parallel algoritlims for the last two problems.
-5-
3. Parallel Approximation of the Jacobian
The approach used in the parallel approximation of the Jacobian by forward differences
depends on two criteria; (1) whether the function F (x) is sufficiently separable, and (2) if the
function is not separable, whether evaluation of the function is computationally expensive.
These criteria are somewhat subjective; whether they apply is a symptom of the specific prob-
lem considered. In this section we describe parallel Jacobian approximation schemes which
deal with the cases delineated by these criteria.
As before, let I? be the set of row indices assigned to processor i and let (x) be this
set of the rows of the Jacobian evaluated at the point x. When evaluation of the function is
separable, then the Jacobian can be evaluated in parallel by having each processor compute its
components of the I -th column according to the formula
Jj?(x)ei = FJj(X+T?j)--H FJj(X)			(3.1)
T
It is often the case that evaluation of F (x) is not completely separable, there may be some
amount of redundant computation due to common factors that must be computed for each parti-
tion of the function Flj(X) i =1p. If this redundant computation is inexpensive relative to
communication cost entailed by using a column-oriented scheme then we consider this compu-
tational overhead tolerable. All of the test problems considered in the experimental section fall
into this category. Otherwise, if the redundant computation required by such a partition of the
rows is deemed too expensive, a column-oriented approach to approximating J(x) must be
adopted. In this approach a set of column indices, Kj, is assigned to each processor i. For
each k E Kj the k -th column of J is approximated at processor i by the usual forward
differences formula,
J(x)4 =--H F(x+Te?)--HF(x)			(3.2)
T
where in this case the function F can be thought of as a "black box." Assuming that the com-
putation required to evaluate F (x) is independent of x, the computational load can be balanced
by making the index sets Kj as close to the same size as possible. Using a wrap mapping will
keep the cardinality of these sets within one, therefore any disparity in workload becomes rela-
tively better as n Ip increases.
-6-
The resulting problem is how to get this column-oriented distribution of the data con-
verted to a row-oriented distribution for the QR factorization stage. Foitunately, there exist
efficient algorithins for transposing a matrix on the hypercube [MvV87]. Of course, for n
large enough this whole problem can be avoided by using column-oriented algorithms in the
first place. Experimentally, we did not take such an approach, but good column-oriented QR
factorization algorithins exist [M87]. In addition, in section 5 we describe an efficient
column-oriented version of the Levenberg-Marquardt algorithins.
A more subtle problem occurs when nip is small, the evaluation of F(x) is not separable,
and the this evaluation is computationally expensive relative to the QR factorization. Suppose
a step Pk is to be considered at the k -th iteration of the algorithin then F (Xk+Pk) must be
evaluated to deterinine if it meets certain acceptance criteria. When this computation is rela-
tively expensive and not separable, and therefore must be done on one processor, then the
remamder of the processors will remain idle during this computation. This can result in detri-
mental effects on the efficiency of the ernire implementation. Byrd, Scnnabel, and Shultz
[BSS88] and Coleman and Li [cL87J note that this problem can be alleviated somewhat by
guessing, based on the previous iteration, whether the proposed point will be accepted. If
acceptance is assuuied, then the Jacobian at Xk+Pk can begin to be approximated by idle proces-
sors. If we guess that the proposed iterate will not be accepted, then idle processors could
evaluate the ftmction at some additional points which might fare better with the acceptance cri-
teria. These ideas were not implemented in our code but could easily be added for cases of
difficult fanctions. However, for nip ? 1, the computation required to estimate the Jacobian
will always dominate these isolated function evaluations.
-7-
4. A New Parallel Row-Oriented Householder QR Algorithm
The efficiency of the parallel QR factorization used to solve equation (2.4) is of
paramount importance because a completely new approximation to the Jacobian is computed
for each iteration. Consequently, a full QR factorization is also required. As we will see in
Section 6 the QR factorization represents the dominant computational cost for the test problems
we considered. In this section we present a new parallel row-oriented Householder QR factori-
zation that was found to be more efficient than previous hybrid (Householder/Givens) factoriza-
tion algorithms. In addition, this algorithin has the advantage that it produces the same House-
holder vectors that would be produced by a standard sequential Householder QR algoritllm
(unlike the hybrid scheme) which can be advantageous in situations where the same system
must be solved for multiple right hand sides. Finally, we show that column pivoting can be
introduced into the algorithin with only a slight increase in the computation and communication
complexity. In our implementation column pivoting is impoflant because the QR factorization
can then be used to estimate matrix rank.
Most of the research on QR algorithins for the hypercube has been directed toward
column-oriented methods, however two row-oriented algorithins have been considered previ-
ously [CP86, PR87]. These two algorithins are very similar in that to reduce each column of
the matrix first a reduction involving only data local to each processor is performed and then
this stage is followed by a global reduction requiring communication between the processors.
The reduction of rows local to a processor results in one row per processor with a norizero in
the column being reduced to upper triangular form. The advantage of this approach is that all
these reductions and matrix updates will be local to the processors and with the wrap mapping
of rows the computational load will be well-balanced. Following this local stage is a global
stage where a minimum spanning tree is embedded in the hypercube rooted at the processor
where the nonzero for the column under consideration should reside. Rows are communicated
up this tree and the leading nonzero is annihilated by a Givens rotation with respect to the
parent's row. These rows are then updated with this rotation and the result communicated
back to the child. The hypercube topology allows this global reduction process to take place in
log (p) steps. Of these two algorithms the one presented by Pothen and Raghavan [PR87j
seems to be the most efficient since Householder reductions, as opposed to Givens, are used in
the local stage.
-8-
Our algorithm is computationally more efficient than the hybrid approach because the full
Householder vector is calculated and the intermediate Givens reductions are avoided. The
difficulty is obtaining the same communication complexity as the hybrid approach. The trick is
to notice that computation of the Householder vector and the rank one update to the matrix can
be combined to half the number of messages that seem to be required at first inspection.
To illustrate the algorithm consider the QR factorization of an m x n matrix A. At step I
of the factorization the first j--Hl rows of R and the Householder vectors have been computed
and we need only consider the (m--Hj+l)x(n--Hj+l) lower right submatrix of A, which we will
denote by A ?), with columns a???, k = jn. The Householder transformation, p(i), to reduce
the first column of A??, a??, is given by
T
p?)?[I?2 V?)V(j) ?
V(j)V(j)
(4.1)
where ?&) = a?? + IIa.&)112 e?. To determine a??+1?, k=j+ln, we need to compute the
corresponding rank one update to A ?) given by:
T2 vT??a?(i) Vcj)			(4.2)
Y(j)V(j)
(4.3)
with ak&) defined as shown. Suppose row j is assigned to a processor that we will designate
leader. Note that vc,? agrees with a ?? except in the first component, hence the portions of the
inner product vT&?a?(i) local to each processor are just a?? Ta?Ci) except on leader where a?? and
V(j) differ in the first component. We can take advantage of this fact to combine the communi-
cation to compute V(i) with the communication required for the rank one update to the
remamder of the matrix. An oufline of the resulting algorithin is given as Algorithin 4.1. For
this description we use the notation [a???i1 to represent the subvector of a?&) with components
given by the index set j?. The a vector is a work vector used in the computation of IIa???II2 and
the constants ?k??, k =j+1
-9
Index Set: Ij (set of rows indices assigned to processor i 1
Proc (i) : (prognim for processor i)
Forj = I,			, n do
If (i =leader) Delete (jJ Irom j?
Compute dot products for k = ----- - n
= [a???J [a???],
Combine [O?j,' , ??] using gather-sum
If (i leader) then
Compute first component of vc'? and the coefficients
- ?)j and broadcast the result
Update columns, k = ....... , n
[a??+1?i1 =			--H ak&) [y?)J1
enddo
Algorithm 4.1: A Panillel Row-Oriented Householder QR Algorithm
.8
.6
? 4
.2
0
Algorithm 4.1
(Householder) ?
Hybrid Algorithm
0			250
500			?750			1000
Number of Rows (m)
1250
Figure 4.1: Efficiency of AIgorithm 4,1 and Hybrid on a 16 Node Hypercube (n=lOO)
- 10 -
In Figure 4.1 we show the efficiencies of this algoritnm compared to the hybrid algoritlirn
described by Pothen and Raghavan in [PR87] as a function of the number of rows. For this
plot the number of columns is fixed at 100. All of the results presented in this paper were
obtained on a 16 node iPSC/l hypercube with extended memory. The efficiencies shown were
calculated by dividing the time taken by an efficient sequential implementation of the algorithiii
run on one processor by the product of the time taken by the parallel implementation and the
number of processors used. In this case our parallel implementations were compared with the
MINPACK QR factorization subroutine QRFAC executed on a single processor of the hyper-
cube and the efficiencies shown were computed from the execution times of these programs.
In Table 4.1 we show some representative execution times for our impl?i'i?iitations of the
hybrid algorithin (Hybrid) and Algoritlim 4.1 (Householder) as compared to the sequential QR
factorization prograni (Single Processor).
Execution Times (sec)			Efficiencies
100			500			434.2			35.7			31.5			.76			.86
Table 4.1: Execution Times and Efficiencies of Algorithin 4.1 on a 16 Node Hypercube
A subtle point in solving equation (2.5) is that the orthogonal matrix Q does not need to
be saved if the right hand side of equation is updated along with the rows of the Jacobian. To
achieve this aim the right hand side is treated like an additional column to the matrix J and
distributed across the processors in the same wrap mapping and updated along with the
corresponding rows of the Jacobian.
- 11 -
Column pivoting can be added to Algorithm 4,1 in the following manner. The column
norms of the matrix A are initialized at the beginning of the algoriflirn and are then updated
after each stage of the computation to obtain the column norms of A (`). For example, suppose
at stage ? in the QR factorization the column norms a?&) 112, k = 1n are known by leader.
The column of maximum norm, krn??, is determined by leader, the result is broadcast, and
columns j and kmax are interchanged by all processors. After stage I of the algorithm the
updated norms can be obtained by the formula:
II ak?+'? 112 = lla?&)ll? 1 --H
lla???ll?
(4.4)
for k =j+l, , n, and the results sent to the next leader (i.e. the next processor on the ring)
for stage j+1 of the QR algorithm. Note that numerical cancellation can be a problem in com-
puting these norm updates. However, circumstances that would result in this problem can be
monitored for and in such cases the suspect column norm can be recomputed. A standard way
to monitor for numerical cancellation is to keep track of the products of the multiplicative fac-
tors m equation (4A) that have been obtained since the last explicit calculation of the column
norm. When this product is sufficiently less than one, then there is the possibility of cancella-
tion error and the column norm is recomputed. Iii uur iiiipl?ni?iitation the recompu?ation is
done by broadcasting a special notifier to the other processors instead of the column pivot.
The required column norms are then recomputed and the result gathered at leader. Our obser-
vation has been that recomputation of the column norms is rarely required and therefore does
not significantly effect the efficiency of the algorithm.
One last concern for numerical stability might be the possibility of overflow from the
way the ?Lk are computed in Algorithm 4.1. We note that these partial sums can be scaled by
the most recent approximation to the column norms available to all the processors. However,
we did not find it necessary to include this scaling in our implementation.
A graph comparing the efficiencies of AIgorithm 4.1 with and without pivoting is
included as Figure 4.2. The efficiency is again computed by comparing the running times of
parallel algorithms to running times of the MINPACK QR subroutine QRFAC on a single pro-
cessor. These results were obtained on a 16 node hypercube with the number columns fixed at
lOO. In table 4.2 we include some representative times from these experiments.
- 12 -
.8
.6
.4
.2
0
-H
-H
Algorithm 4.1
without pivoting
Algorithm 4.1
with pivoting
0			250
500			?0			1000
Number of Rows (m)
1250
Figure 4?: Efficiency of Algoritllm 4.1 With and Without Column Pivoting
Execution Times (sec)			Efficiencies
100			500			439.0			34.1			.80
100			1000			906.5			62.1			.91
Table 4.2: Execution Times and Efficiencies of Algorithin 4.1 with Pivoting
- 13 -
5. A Parallel Implementation of the Levenherg-Marquardt Algorithm
To determine the Levenberg-Marquardt parameter Algorithin 5.1 was used to reduce
equation (2.7) to upper triangular form. This reduction requires a surprising amount of work:
n (n+l)/2 Givens rotations and the corresponding row updates, or 0(n3) flops. Note that the
work required in this reduction is independent of m, the number of rows. AIgorithm 5.1
details a parallel method to accomplish this reduction. In the algorithmic description let S
represent storage for an upper triangular matrix which is initially set equal to the matrix ?I in
equation (2.7). Remember that the rows of R and S are wrapped onto an embedded ring of
processors as described in Section 1.
Functions: next (returns number of next processor in the ring),
prev (returns number of previous processor in the ring)
Index Set: !? (set of rows assigned to processor i)
Proc (i) : (program for processor i)
Forj =O,n--Hl do
If (j?O) receive rows skTH, k e1), from processor prev(i)
Fork?i?? do
Compute Clivens rotation to zero the bottom of the vector (Rk? Sk?j .k
Update rows RkT and SATH with above Givens rotation
If (k=]+l) Delete (k) from I?
enddo
Send rows SkT-j, k E i???, to processor next(i)
enddo
Algorithm 5.1: A Parallel Row-Oriented R-S Reduction
Algorithm 5.1 proceeds in n stages which have been indexed by 1=0n--Hl in the
description. At stage j of algorithm 5.1 the superdiagonal of S that is a distance j from the
main diagonal is eliminated by Givens rotations. After n stages, the upper triangular matrix S
has been completely zeroed and the updated upper triangular matrix R is still wrapped onto the
processors in the same manner as at the stan of the algoritlim. As the leading nonzero of each
row of S is eliminated and the corresponding rows updated the rows of S move around the
embedded ring in a systolic manner. Note, that although the work at each stage is not com-
pletely balanced, the processor doing the most work rotates around the ring. This imbalance is
14
somewhat offset by the required communication. Experimental results of the efficiency of this
algorithm as a function of number of columns are shown in Figure 5.1.
.8
.6
.4
.2
0
(row--Horiented)
Algorithm 5.2
(column-oriented)
0			100
200			300			400
Number of Columns (n)
Figure 5.1: Efficiencies of Algorithms 5.1 and 5.2 on a 16 Node Hypercube
500
A column-oriented approach is also possible and is presented as Algorithm 5.2. From
Figure 5.1 it is apparent that this version is superior experimentally. Even though the algo-
rithm is a bit more complicated, the total number of messages that have to be sent is the same
as in Algorithin 5.1 and the total number of values that have to communicated is actually less.
For an average step j in Algorithin 5.2 we need only communicate the single Givens vector g
of length 0(n) between neighboring processors. For the row-oriented version we need to com-
municate 0 (nip) rows of S of length 0(n) between processors. In practice the rows of S are
combined into one long message which results the same number of communication start-ups as
appear in Algorithin 5.2. The message start-up cost, measured in "flops," for the Intel iPSC
hypercube is very expensive and is normally the dominant factor in the communication cost of
an algorithm. However, for these two algorithms the average message lengths are extremely
different, hence one sees a difference in efficiency between these two algorithms experimen-
tally.
- 15 -
Index Set: Kj (set of columns assigned to processor i)
Proc (i) : (program for processor i)
Forl =O,n--Hl do
If (j?<)) Receive Clivens vector g from processor prev(i)
F0rkEKj do
Fort =min(j,p--Hl),..,ldo
Update rows RkT?I and SkT?j with Givens rotation g?-?
enddo
Compute (livens rotation to zero the bottom of the vector (Rk.k , Sk?j?)T
Update rows RAT and SkT?j with above (livens rotation
Update g? in Givens vector
If (k=j+l) Delete (k) from I?
enddo
Send Givens vector g to processor next(i)
enddo
Algorithm 5.2: A Parallel Column-Oiiented R-S Reduction
For Algoritlirn 5.2 the columns, as opposed to the rows, of R and S are wrapped onto the
ring of processors. Rather than communicating rows of S between neighboring processors, the
Givens rotations are stored in vectors g which rotate around the ring. Once the algorithm has
been running for more than p steps, i.e. j >--Hp--Hl, then the Givens vector g is completely filled
with updates that need to be applied once received. The order that these rotations are applied
in the I loop is important. Since they operate on the same row of S, the rotations must be
applied from oldest to newest. Also, note by row RAT we mean the nonz:ero components of row
RAT that are local to processor i. These components are given by the index set Kj. Note that
the difference in the average length of the messages sent between Algoritlims 5.1 and 5.2
translates into a much more efficient algorithni for the column-oriented approach. It is possible
to decrease the length of the messages sent in the row-oriented method by postponing the
application of the Givens rotations. Unfortunately, this algorithm is very complicated and was
not implemented.
The reduction of equation (2.7) to upper triangular form is the major task to be rendered
parallel in an algoritlim for determining the Levenberg-Marquardt parameter X* and we have
- 16
shown that there exist effective algorithins to perform this reduction. However, efficient solu-
tion of triangular systems is also important in this context. In fact, for each iteration involving
a solution of equation (2.7) there are two associated triangle solves that are used to bracket the
solution ? [M78]. Experirnentally we used the triangle solve algorithms developed by Li and
Coleman [LC87, LC88], but it should be noted that the efficiencies of these algorithms are not
neafly as good as those of Algoritlirns 5.1 and 5.2. This difference is what accounts for the
discrepancies between the efficiencies shown in Figure 5.1 and the efficiences reported in the
next section for solving for the Levenberg-Marquardt parameter. This effect is apparent even
though there is an 0(n) difference between the amount of work required for AIgorithm 5.1 and
the corresponding triangular solves. The importance of efficient parallel triangle solvers has
also been observed in the solution of systems of nonlinear equations [CL87].
- 17 -
6. Experimental Results and Conclusions
These algoritlirns were implemented on a 16-node Intel iPSC/l hypercube with extended
memory in RMfFoflran and run under version R3.l.l of the iPSC operating system. The
efficiencies shown below were calculated by dividing the running time of MINPACK
[MGH8O] code on a single processor by 16 times the running times of the algorithins
described above on the hypercube. The comparison is fair as both programs generate the same
sequence of iterates and consequently do the same number of Jacobian approximations, QR
factorizations, and Newton iterations in computing the Levenberg-Marquardt parameten
Problem			F Eval.			J Approx.			F Eval.
Function Form			Characteristics
Number			Cost			Cost			Separable?
A E R???, full rank,
1			F A x --He			A x easily evaluated,			0(m)			cheap			no
e an rn-vector of ones
A E Rmxn, low rank,
2			F=Ax--He			0 (mn)			expensive			yes
dense
3			fi =?Aq?XjX??l			A sparse			0(m)			cheap			yes
j,k
fi =?(CXP(xjTi)			Constrained, x <0 for
4			j			0 (nut)			expensive			yes
?eXP((?jTi))			(x<O and T>O
Table 6.1: Description of Test Functions
The test problems used to obtain the experimental results are described in Table 6.1.
Shown is the functional form of the test problems and the computational complexity of evaluat-
ing each function is given in the column labeled "F Eval. Cost." We also make a subjective
determination as to whether estimation of the Jacobian is cheap or expensive relative to its QR
factorization. To be more specific about the functions used for testing, problem 1 has an 0(m)
evaluation cost because of a special form for the matrix A: Ajj =--H2/m, i ?j; Ajj = 1 --H2/m. The
matrix A used in problem 2 is of low rank: Aq = ii +0 (? perturbations, where e is the machine
- 18 -
precision. For problem 3 the tensor A has a bandwidth of 5 allowing for flinction evaluation
with 0(m) cost The constants used for problem 4 were (xj = --Hj and Tj = jim In the final
column of Table 6.1 we note whether we consider evaluation of these flinctions to be separ-
able.
Efficiencies Compared to MINPACK
(% time spent in routine)
Prob			n			m			Total			QR			L-M			R-S			Tri. S.			J Appr
.644			.654			---			---			.089			.322
1			100			250
(100.0)			(95.8)			(0.0)			(0.0)			(1.2)			(2.8)
.801			.823			---			---			.089			.322
1			100			500
(100.0)			(94.8)			(0.0)			(0.0)			(1.4)			(3.8)
.679			.649			.527			.631			.084			.821
2			100			250
(100.0)			(46.2)			(23.7)			(19.1)			(4.6)			(28.8)
.784			.821			.510			.634			.093			.822
2			100			500
(100.0)			(48.1)			(15.2)			(11.8)			(2.8)			(35.6)
.813			.823			.758			.806			.175			.817
2			200			500
(100.0)			(44.5)			(19.2)			(17.7)			(1.2)			(35.7)
.798			.824			.547			.626			.098			.937
3			100			500
(100.0)			(70.5)			(14.2)			(12.0)			(2.2)			(11.4)
.814			.826			.759			.807			.175			.940
3			200			500
(100.0)			(50.6)			(39.6)			(36.6)			(2.5)			(9.1)
.824			.823			.756			.810			.179			.859
4			200			500
(100.0)			(36.1)			(25.0)			(23.1)			(1.6)			(38.5)
Table 6.2: Experimental Results of Parallel Algorithins Compared with MINPACK
- 19 -
Table 6.2 summerizes the experimental results obtained by comparing our parallel algo-
rithlns with the MINPACK code running on a single processor for solving the test problems
described in Table 6.1. The efficiences and the fraction of the total parallel running time spent
in each of six sections of the programs are detailed. The six sections refer to: QR, the QR
factorization of the Jacobian approximation; L-M, computation of the Levenberg-Marquardt
parameter; R-S, the R-S reduction described in Section 5; Tri. 5., the solution of triangular sys-
tems; and J Appr., the approximation of the Jacobian by forward differences.
These results were choosen to illustrate an average case behavior of the parallel algo-
rithins in solving nonlinear least-squares problems. For a particular problem the fraction of
time spent in the different routines and the number of iterations required for convergence can
vary dramatically for various choices of a starting point, initial trust-region size, and termina-
tion tolerances. For the above problems the MINPACK tolerances were set to ?, where ? is
the machine precision. The initial trust-region size was given by 10011 io 112 where i0 is the
starting point normalized by the 2-norms of the columns of the initial Jacobian approximation.
For problems 1, 2, and 3 an n-vector of ones was used as the starting point, and for problem 4
the starting point Xj = --H (I + 0.1) was employed. For the problem sizes shown in Table 6.2
problem 1 required respectively 3 and 4 iterations (Jacobian approximations) for convergence,
problem 2 used 2, 3, and 3 iterations, and problem 3 required 6 and 14 iterations to converge.
More than 50 iterations were required for problem 4 to converge; the results shown in Table
6.2 were taken from the first 10 iterations.
From these results it is apparent that the QR factorization is the dominant computational
cost for these problems. Consequently, the implementation is more efficient as m, the number
of rows, increases. The cost of computing the Levenberg-Marquardt parameter can also be
significant; this computation could dominate the QR factorization for problems that require a
disproportionately large number of R-S reductions and for which the ratio min is close to one.
This effect would occur in problems where the Jacobian is rank deficient or the function is
very nonlinear (which would require a small trust-region in order to ensure an accurate qua-
dratic model of the function). But again, for a fixed ratio min, the efficiency in solving for the
Levenberg-Marquardt parameter increases as m increases. As expected the parallel triangular
system solves are very inefficient when compared to the QR factorization and R-S reductions.
However, the time required for these solves comprises only a small fraction of the total
- 20 -
computation time of the parallel implementation and therefore it does not significantly decrease
the overall efficiency.
For functions whose evaluation is very expensive we observe that the computation
required for the approximation of the Jacobian by forward differences can equal or exceed the
computation required for these other tasks. For the expensive test functions we considered the
function evaluation was separable and hence the row-oriented Jacobian approximation algo-
rithm yielded efficiencies comparable to the QR factorization and Levenherg-Marquardt param-
eter solves. If the function evaluation were expensive and not separable one would have to
resort to a column-oriented Jacobian approximation algorithm as described in Section 3. In
this case the efficiency of the implementation would depend on the number of columns heing
sufficiently large. Another possibility would be to obtain an algebraic expression for the Jaco-
bian and evaluate the Jacobian directly rather than using forward differences; the algebraic
evaluation of the Jacobian could be separable even when the function itself is noL
In sum, we have observed good efficiencies when solving moderately sized nonlinear
least-squares problems on the Intel hypercube. For the test problems considered we noted that
the efficiency of our parallel implementation inipi?v?d as the ratio mip increased where m is
the number of rows of Jacobian and p is the number of processors. We also point out that it
is possible to solve much larger problems than those we have described above (which had to
be run on one processor for comparison). The efficiencies for such larger problems would be
correspondingly better. An additional topic for future consideration would be the solution of
large sparse problems, especially problems that have special structure. For example, we note
that a simple generalization of Algorithm 5.1 would work well if the Jacobian were banded and
that the row-oriented approach to the QR factorization would be efficient in solving problems
with a block structure.
Acknowledgements
The work that produced this paper was partially completed with the assistance of comput-
ing facilities of the Advanced Computing Facility at the Cornell Center for Theory and Simula-
tion in Science and Engineering, which is supported by the National Science Foundation and
New York State. We would also like to acknowledge discussions at ORNL (as part of their
Numerical Linear Algebra Year).
- 21
References
[BSS88] R. Byrd, R. Schnabel, and G. Shultz. Parallel Quasi-Newton Methods for Uncon-
strained Optimization. Department of Computer Science, University of Colorado
at Boulder, 1988.
[CP86] Richard Chamherlain and M. J. D. Powell. QR Factorisation for Linear Least-
Squares Problems on the Hypercube. Chr. Michelsen Institute, 1986.
[C84J Thomas F. Coleman. Large Sparse Numerical Optimization. Springer-Verlag,
1984.
[CL87] Thomas F. Coleman and Guangye Li. Solving Systems of Nonlinear Equations on
a Message-Passing Multiprocessor. Tech. Rep. 87-887, Computer Science Depart-
ment, Cornell University, 1987.
[LC87] Guangye Li and Thomas F. Coleman. A New Method for Solving Triangular Sys-
tems on Distributed Memory Message-Passing Multiprocessors. TR 87-812, Com-
puter Science Department, Cornell University, 1987.
[LC88j Guangye Li and Thomas F. Coleman. A Parallel Triangular Solver for a Distri-
buted Memory Multiprocessor. SlAM J. Sci. Comput., 1988.
[M87j C. Molen Matrix Computations on Distributed Memory Multiprocessors. Tech.
Rep., Intel Scientific Computers, 1987.
[MVV87] 0. M. McBryan and E. F. Van de Velde. Hypercube Algorithms and Implementa-
tions. SlAM J. Sci. Comput., 1987.
[GMW8l] Phillip Gill, Walter Murray, and Margaret Wright. Practical Optimization.
Academic Press, 1981.
[M77] Jorge J. Mor?. The Levenberg-Marquardt Algorithm: Implementation and Theory.
Springer-Verlag, 1977.
[MGH8O] Jorge J. Mor?, Burton Garbow, and Kenneth Hillstrom. User Guide for
MINPACK-1. Argonne National Laboratory, 1980.
[PR87] Alex Pothen and Padma Raghavan. Distributed Orthogonal Factorization.' Givens
and Householder Algorithms. Department of Computer Science, The Pennsylvania
State University, 1987.