Chapter 6
NON-LINEAR REGRESSION MODELS
TAKESHI AMEMIYA*
Stanford University
Contents
1. Introduction
334
2. Single equation-i.i.d. case
336
2.1. Model
336
2.2. Asymptotic properties
337
2.3. Computation
341
2.4. Tests of hypotheses
347
2.5. Confidence regions
352
3.
Single equation-non-i.i.d. case
354
3.
I.
Autocorrelated errors
354
3.2. Heteroscedastic errors
358
4. Multivariate models
359
5. Simultaneous equations models

econometrician to estimate an increasing number of non-linear regression models
in recent years. Non-linearity arises in many diverse ways in econometric applica-
tions. Perhaps the simplest and best known case of non-linearity in econometrics
is that which arises as the observed variables in a linear regression model are
transformed to take account of the first-order autoregression of the error terms.
Another well-known case is the distributed-lag model in which the coefficients on
the lagged exogenous variables are specified to decrease with lags in a certain
non-linear fashion, such as geometrically declining coefficients. In both of these
cases, non-linearity appears only in parameters but not in variables.
More general non-linear models are used in the estimation of production
functions and demand functions. Even a simple Cobb-Douglas production
function cannot be transformed into linearity if the error term is added rather
than multiplied [see Bodkin and Klein (1967)]. CES [Arrow, Chenery, Minhas and
Solow (196 l)] and VES [Revankar (1971)] production functions are more highly
non-linear. In the estimation of expenditure functions, a number of highly
non-linear functions have been proposed (some of these are used in the supply
side as well)-Translog [Christensen, Jorgenson and Lau (1975)], Generalized
Leontief [Diewert (1974)], S-Branch [Brown and Heien (1972)], and Quadratic
[Howe, Pollack and Wales (1979)], to name a few. Some of these and other papers
with applications will be mentioned in various relevant parts of this chapter.
The non-linear regression models I will consider in this chapter can be written
in their most general form as
(1.1)
where y,, .x,, and (Y~ are vectors of endogenous variables, exogenous variables, and
parameters, respectively, and uif are unobservable error terms with zero mean.
Eqs. (1. l), with all generality, constitute the non-linear simultaneous equations
model, which is analyzed in Section 5. I devote most of the discussion in the
chapter to this section because this area has been only recently developed and
therefore there is little account of it in general references.
Ch. 6: Non -linear Regression Models

in this Handbook). There are a few other important topics which, although
non-linearity is involved, woud best be studied within another context, e.g.
non-linear error-in-variable models and non-linear time-series models. Regarding
these two topics, I recommend Wolter and Fuller (1978) and Priestley (1978).
Finally, I conclude this introduction by citing general references on non-linear
regression models. Malinvaud (1970b) devotes one long chapter to non-linear
regression models in which he discusses the asymptotic properties of the non-
linear least squares estimator in a multivariate model. There are three references
which are especially good in the discussion of computation algorithms, confidence
regions, and worked out examples: Draper and Smith (1966) Bard (1974) and
Judge, Griffiths, Hill and Lee (1980). Several chapters in Goldfeld and Quandt
(1972) are devoted to the discussion of non-linear regression models. Their
Chapter 1 presents an excellent review of optimization techniques which can be
used in the computation of both the non-linear least squares and the maximum
likelihood estimators. Chapter 2 discusses the construction of confidence regions
336
T. Amemiya
in the non-linear regression model and the asymptotic properties of the maximum
likelihood estimator (but not of the non-linear least squares estimator). Chapter 5
considers the Cobb-Douglas production function with both multiplicative and
additive errors, and Chapter 8 considers non-linear (only in variables) simulta-
neous equations models. There are two noteworthy survey articles: Gallant
(1975a), with emphasis on testing and computation, and Bunke, Henscheke,
Strtiby and Wisotzki (1977), which is more theoretically oriented. None of the
above-mentioned references, however, discusses the estimation of simultaneous
equations models non-linear both in variables and parameters.
2.
Single equation-i.i.d. case
2.1.
Model

The non-linear least squares (NLLS) estimator, denoted p, is defined as the
value of /I that minimizes the sum of squared residuals
S,(P) = t
[Yt
-fhP)12.
(2.5)
It is important to distinguish between the p that appears in (2.5), which is the
argument of the function f(x,, m), and &, which is a fixed true value. In what
follows, I will discuss the properties of p, the method of computation, and
statistical inference based on 8.
2.2.
Asymptotic properties
2.2.1. Consistency
The consistency of the NLLS estimator is rigorously proved in Jennrich (1969)
and Malinvaud (1970a). The former proves strong consistency (j? converging to
&, almost surely) and the latter weak consistency (p converging to &, in
probability). Weak consistency is more common in the econometric literature and
is often called by the simpler name of consistency. The main reason why strong
consistency, rather than weak consistency, is proved is that the former implies the
latter and is often easier to prove. I will mainly follow Jennrich’s proof but
translate his result into weak consistency.
The consistency of b is proved by proving that plim T- ‘S,( j3) is minimized at
the true value &. Strong consistency is proved by showing the same holds for the
almost sure limit of
T- ‘S,( /3)
instead. This method of proof can be used to prove
the consistency of any other type of estimator which is obtained by either
minimizing or maximizing a random function over the parameter space. For
example, I used the same method to prove the strong consistency of the maximum
likelihood estimator (MLE) of the Tobit model in Amemiya (1973b).

P, PUT-‘&-(P)- S(P)1 ’ &I< 6.
It is easy to construct examples in which the violation of any single assumption
above leads to the inconsistency of 8. [See Amemiya (1980).]
I will now give a sketch of the proof of the consistency and indicate what
additional assumptions are needed as I go along. From (2.1) and (2.5), we get
=A,+A,+A,,
(2.9)
where c means CT=, unless otherwise noted. First, plim A, = ut by a law of large
numbers [see, for example, Kolmogorov Theorem 2, p. 115, in Rao (1973)].
Secondly, for fixed &, and p, plim A, = 0 follows from the convergence of
T-‘C[f,(&)- f,(p)]’
by Chebyshev’s inequality:
Since the uniform convergence of A, follows from the uniform convergence of the
right-hand side of (2.10), it suffices to assume
converges uniformly in
fi, ,
& E
B.
(2.11)
Having thus disposed of A, and
A,,
we
need only to assume that lim
A,
is
uniquely minimized at PO; namely,
lim+E[f,(&)-N)l’-o ifP*&.
(2.12)
To sum up, the non-linear least squares estimator B of the model (2.1) is
consistent if (2.6), (2.1 l), and (2112) are satisfied. I will comment on the signifi-

exists and is positive definite. It can be easily proved that
in the linear model the above assumption is not necessary for the consistency of
least squares and it is sufficient to assume
(X’X)-
’
+ 0. This observation
suggests that assumption (2.12) can be relaxed in an analogous way. One such
result can be found in Wu (198 1).
2.2.2.
Asymptotic normality
The asymptotic normality of the NLLS estimator B is rigorously proved in
Jennrich (1969). Again, I will give a sketch of the proof, explaining the required
assumptions as I go along, rather than reproducing Jennrich’s result in a theo-
rem-proof format.
The asymptotic normality of the NLLS estimator, as in the case of the MLE,
can be derived from the following Taylor expansion:
(2.13)
where a2$/apap’ is a K
x
K
matrix of second-order derivatives and p* lies
between j? and &. To be able to write down (2.13), we must assume that f, is
twice continuously differentiable with respect to p. Since the left-hand side of
(2.13) is zero (because B minimizes S,), from (2.13) we obtain:
@(~_p,)=_
1
a2sT
[
Twl,.]‘$ %I,,-
(2.14)

we assume
exists and is non-singular,
then
1
as
t
0
afi
PO
+ N(0,4&).
(2.17)
(2.18)
This result can be straightforwardly obtained from the Lindberg-Feller central
limit theorem [Rao (1973, p. 128)] or, more directly, from of Anderson (197 1,
Theorem 2.6.1, p. 23).
Proving (ii) poses a more difficult problem. Write an element of the matrix
~-l(a~s,/apap)~.
ash@*).
0
ne might think that plim hT( /3*) = plim hT( &,)
follows from the well-known theorem which says that the probability limit of a
continuous function is the function of the probability limit, but the theorem does
not apply because h, is in general a function of an increasing number of random
variables y,, j2,.
. . ,y,.
But, by a slight modification of lemma 4, p. 1003, of
Amemiya (1973b), we can show that if hr( p)
converges almost surely to a certain
non-stochastic function h( /?) uniformly in p, then plim hT( p*) = h(plim /I*) =
h( &). Differentiating (2.15) again with respect to p and dividing by

a F
x K
matrix. Note that (2.24) exactly holds
in the linear case. The practical consequence of the approximation (2.24) is that
all the results for the linear regression model are asymptotically valid for the
non-linear regression model if we treat G as the regressor matrix. In particular, we
can use the usual t and
F
statistics with an approximate precision, as I will
explain more fully in Sections 2.4 and 2.5 below. Since the matrix G depends on
the unknown parameters, we must in practice evaluate it at b.
2.3.
Computation
Since there is in general no explicit formula for the NLLS estimator b, the
minimization of (2.5) must usually be carried out by some iterative method. There
342
T. Amemiya
are two general types of iteration methods: general optimization methods applied
to the non-linear least squares problem in particular, and procedures which are
specifically designed to cope with the present problem. In this chapter I will
discuss two representative methods - the Newton-Raphson iteration which be-
longs to the first type and the Gauss-Newton iteration which belongs to the
second type - and a few major variants of each method. These cover a majority of
the iterative methods currently used in econometric applications. Although not
discussed here, I should mention another method sometimes used in econometric
applications, namely the so-called conjugate gradient method of Powell (1964)
which does not require the calculation of derivatives and is based on a different
principle from the Newton methods. Much more detailed discussion of these and
other methods can be found in Chapter 12 of this Handbook and in Goldfeld and
Quandt (1972, ch. 1).

where
I
is the identity matrix and (Y, is a scalar to be appropriately chosen by the
researcher subject to the condition that ( a2&/apap’)jn + a,Z is positive definite.
This modification was proposed by Goldfeld, Quandt and Trotter (1966) and is
called
quadratic hill-climbing
(since they were considering maximization). See the
same article or Goldfeld and Quandt (1972, ch. 1) for a discussion of how to
choose (Y, and the convergence properties of the method.
The second weakness may be remedied by the modification:
(2.29)
where the scalar X, is to be appropriately determined. See Fletcher and Powell
(1963) for a method to determine h, by a cubic interpolation of S,(p) along the
current search direction. [This method is called the DFP iteration since Fletcher
and Powell refined the method originally proposed by Davidon (1959).] Also, see
Berndt, Hall, Hall and Hausman (1974) for another method to choose A,.
Ordinarily, the iteration (2.26) is to be repeated until convergence takes place.
However, if B, is a consistent estimator of & such that @(b, - &,) has a proper
limit distribution, the second-round estimator 8, has the same asymptotic distri-
bution as B. In this case, a further iteration does not bring any improvement so
far as the asymptotic distribution is concerned. This is shown below.
By a Taylor expansion of (a&/a/3);, around &, we obtain:
as, as,
-I
I
ab j, =
T
Bo +
where p* lies between B,

non-overlap-
ping consecutive subsets !P,, !P2;,
. . . ,
!PK,
each of which contains
m
elements. If we
define
j$, =m-'ClsuyI
and
f&(/3)=m-LC,,yfi(/3),
i=l,2, ,K, the
Harley-Booker estimator is defined as the value of b that satisfies
K
equations:
Y(i)=(i)(P),
i=1,2
K.
,**.,
(2.34)
Since (2.34) cannot generally be solved explicitly for p, one still needs an
iteration to solve it. Hartley and Booker propose the minimization of EYE ,[ jjCi, -
fCi,(p)12 by an iterative method, such as one of the methods being discussed in
this section. This minimization is at least simpler than the original minimization
of (2.5) because the knowledge that the minimand is zero at /3 = & is useful.
However, if there are multiple solutions to (2.34), an iteration may lead to the
wrong solution.
Hartley and Booker proved the consistency of their estimator. Jennrich (1969)
gave a counterexample to their consistency proof; however, their proof can easily
be modified to take account of Jennrich’s counter-example. A more serious

left-hand side is treated as the dependent variable and (af,/J/3’);, as the vector
of independent variables. Eq. (2.38) reminds us of the point raised above: namely,
the non-linear regression model asymptotically behaves like the linear regression
model if we treat (af/&‘fi’)j as the regressor matrix.
346
T A men1iya
The Gauss-Newton iteration suffers from weaknesses similar to those of the
Newton-Raphson iteration: namely, the possibility of a total or near singularity
of the matrix to be inverted in (2.37), and the possibility of too much or too little
change from /$, to &,+ ,.
In order to deal with the first weakness, Marquardt (1963) proposed a modifi-
cation:
where (Y, is a positive scalar to be appropriately determined by a rule based on the
past behavior of the algorithm.
In order to deal with the second weakness, Hartley (1961) proposed the
following modification. First, calculate
(2.40)
and, secondly, choose A, so as to minimize
$4 ii + %A,)~
O_Ih,_Il.
(2.41)
Hartley proves that under general conditions his iteration converges to a sta-
tionary point: that is, a root of the normal equation &S,/ap = 0. He also proves
(not so surprisingly) that if the iteration is started at a point sufficiently close to
b, it converges to b. See Tomheim (1963) for an alternative proof of the
convergence of the Hartley iteration. Some useful comments on Marquardt’s and
Hartley’s algorithms can be found in Gallant (1975a). The methods of determin-
ing A, in the Newton-Raphson iteration (2.29) mentioned above can be also
applied to the determination of A,, in (2.41).
Jennrich (1969) proves that if the Gauss-Newton iteration is started at a point

h, + 0 almost surely for all n as
T + CO.
(2.W
But (2.44) implies two facts. First, the iteration converges to a stationary point,
and secondly, this stationary point must lie sufficiently close to the starting value
8, since
(,i~-p,)l(~~ ,)~S’6(1+h,+hlXz+ * +h,X,.+_*),
(2.45)
where 6 = & - 8,. Therefore, this stationary point must be B if 8, is within a
neighborhood of & and if b is the unique stationary point in the same neighbor-
hood.
In closing this section I will mention several empirical papers in which the
above-mentioned and related iterative methods are used. Bodkin and Klein (1967)
estimated the Cobb-Douglas (2.2) and the CES (2.3) production functions by the
Newton-Raphson method. Charatsis (1971) estimated the CES production func-
tion by a modification of the Gauss-Newton method similar to that of Hartley
(1961) and showed that in 64 samples out of 74, it converged in six iterations.
Mizon (1977), in a paper the major aim of which was to choose among nine
production functions, including the Cobb-Douglas and CES, used the conjugate
gradient method of Powell (1964). Miion’s article is a useful compendium on the
econometric application of various statistical techniques such as sequential test-
ing, Cox’s test of separate families of hypotheses [Cox (1961, 1962)], the Akaike
Information Criterion [Akaike (1973)], the Box-Cox transformation [Box and
Cox (1964)], and comparison of the likelihood ratio, Wald, and Lagrange multi-
plier tests (see the end of Section 2.4 below). Sargent (1978) estimates a rational
expectations model (which gives rise to non-linear constraints among parameters)
by the DFP algorithm mentioned above.
2.4.
Tests of hypotheses
In this section I consider tests of hypotheses on the regression parameters p. It is

we can find a K,
X K
matrix
R
such that
(R’, Q’) = A’
is non-singular. If we
define (Y = A/3 and partition (Y’ = (‘Y;,), (Y{*)),
the hypothesis Q/3 = c is equivalent
to the hypothesis a(Z) = c.
As noted after eq. (2.24), all the results of the linear regression model can be
extended to the non-linear model .by treating G = ( af/&3’),0 as the regressor
matrix if the assumptions of Section 2.2 are satisfied. Since &, is unknown, we
must use G =
(af/ap)j
in practice. We will generalize the
t
and
F
statistics of
the linear model by this principle. If K, = 1, we have approximately
-qK($z,-%J _
t(T_
K)
gm
’
(2.46)
where L! is the last diagonal element (if &) is the i th element of p, the
i
th diagonal

(2.48)
For each of the four parameters, the empirical distribution of the left-hand side of
(2.46) matched the distribution of
t(T - K)
reasonably well, although, as we
would suspect, the performance was the poorest for 8,.
In testing & = p(Z) when
K, 2
1, we may alternatively use the asymptotic
approximation (under the null hypothesis):
(T-K)[%(i+&@)l _
J-(K
K2wv
29
T_ K)
(2.49)
where b is the constrained non-linear least squares estimator obtained by mini-
mizing S,( /3) subject to pc2,
= pc2). Although, as is well known, the statistics (2.47)
and (2.49) are identical in the linear model, they are different in the non-linear
model.
The study of Gallant (1975~) sheds some light on the choice between (2.47) and
(2.49). He obtained the asymptotic distribution of the statistics (2.47) and (2.49)
under the alternative hypothesis as follows. Regarding S,(b), which appears in
both formulae, we have asymptotically:
S,(B) = u’[l-G(G’G)-‘G’]u,
(2.50)
where G = ( af/a/3’)s, as before. Define G, = ( af/a&,),, Then, Gallant shows
(asymptotically) that
s,(p) = (u+ a),

Now
I consider the test of a non-linear hypothesis
h(P) = 0,
(2.53)
where
h
is a q-vector valued non-linear function such that
q < K.
If /3 are the parameters that characterize a concentrated likelihood function
L(p), where
L
may or may not be derived from the normal distribution, we can
test the hypothesis (2.53) using one of the following well-known test statistics: the
likelihood ratio test (LRT), Wald’s test [WaId (1943)], or Rao’s test [Rao (1947)]:
LRT=2[logL(j)-logL@)], (2.54)
and
(2.55)
(2.56)
‘In deriving the asymptotic approximations (2.51) and (2.52), Gallant assumes that the “distance”
between the null and alternative hypotheses is sufficiently small. More precisely, he assumes that there
exists a sequence of hypothesized values @&) and hence a sequence (/36:> such that fi( &)a - p&)
and fl(&, -P;,;)
converge to constant vectors as T goes to infinity.
*Actually, the powers of the two tests calculated either from the approximation or from the
empirical distribution are identical in testing j3, = 0. They differ only in the test of & = - 1.
Ch. 6: Non -lineur Regression Models
351
where B is the unconstrained maximum likelihood estimator and /? is the
constrained maximum likelihood estimator obtained maximizing L(p) subject to
(2.53).3 By a slight modification of the proof of Rao (1973) (a modification is

3See Silvey (1959) for an interpretation of Rao’s test as a test on Lagrange multi
P
Iiers.
41f 6 is distributed as a q-vector N(0, V), then (5 + a)‘V-‘(6 + p) - x*(q,p’V- p).
‘In the following derivation I have omitted some terms whose probability limit is zero in evaluating
@‘(6’log L/a/3’) and T-‘(a* log L/6’ga/F).
352
T. Amemiya
using a proof similar to Rao’s, we can show that the statistics (2.58), (2.59) and
(2.60) are asymptotically distributed as x’(q) even if u are not normal. Thus,
these statistics can be used to test a non&near hypothesis under a
situation.
In the linear regression model we can show Wald 2 LRT 2 Rao
and Savin (1977)]. Although the inequalities do not exactly hold for
ear model, Mizon (1977) found Wald 2 LRT most of the time in his
2.5.
Confidence regions
Confidence regions on the parameter vector p or its subset can be
constructed
using any of the test statistics considered in the preceding section. In this section I
discuss some of these as well as other methods of constructing confidence regions.
A 100
X
(1 -
(u) percent confidence interval on an element of p can be obtained
from (2.46) as
non-normal
[see Bemdt
the non-lin-
samples.

have the same asymptotic distribution- F(
K, T - K).
I have not come across any
reference discussing the comparative merits of the two methods.
(2.62)
(2.63)
Beale (1960) shows that the confidence region based on (2.63) gives an accurate
result - that is, the distribution of the left-hand side of (2.63) is close to F(
K, T -
K)-
if the “non-linearity” of the model is small. He defines a measure of
Ch. 6: Non
-linear
Regression Models
353
non-linearity as
‘j= ii
2
[.htbi)-h(8)-
~l,(bi-/l)]z.K(r-K)-‘s,(iR)
i=l
f=l
’ { ig, [
,$,
[hcPi~~/,(~,l’]‘)‘~
(2.64)
where
b
b
,, 2,. . . , b,,, are

(2.66)
where 2 is an appropriately chosen
T
X K
matrix of constants with rank K. The
computation of (2.66) is more difficult than that of (2.65) because p appears in
both the numerator and denominator of (2.66). In a simple model where f,( /3) =
P, + Pse%
Hartley suggests choosing Z such that its tth row is equal to
(1,
x,,
xf
).
This suggestion may be extended to a general recommendation that we
should choose the column vectors of Z to be those independent variables which
we believe best approximate G. Although the distribution of the left-hand side of
(2.66) is exactly
F(K, T - K)
for any Z under the null hypothesis, its power
depends crucially on the choice of Z.
354
T. Amemiya
3. Single equation-non4.i.d. case
3. I.
Autocorrelated errors
In this section we consider the non-linear regression model (2.1) where {u,} follow
a general stationary process
cc
U, =
C

Since
A,
involves the vector product f’u and since E(f’u)* = f’Zf $ f’fx,(Z), where
h,(E) is the largest characteristic root of E, assumption (2.11) implies plim
A, = 0
by Chebyshev’s inequality, provided that the characteristic roots of 2 are bounded
from above. But this last condition is implied by assumption (3.3).
To prove the asymptotic normality in the present case, we need only prove the
asymptotic normality of (2.16) which, just as in the linear model, follows from
theorem 10.2.11, page 585, of Anderson (1971) if we assume
I5
IYjl
<
O”
(3.4)
j=O
in addition to all the other assumptions. Thus,
~(B-Po)-~[0,0,21imT-‘(G’G)-‘G’~G(G’G)~’],
(3.5)
Ch. 6: Non
-linear
Regression Models
355
which indicates that the linear approximation (2.24) works for the autocorrelated
model as well. Again it is safe to say that all the results of the linear model are
asymptotically valid in the non-linear model. This suggests, for example, that the
Durbin-Watson test will be approximately valid in the non-linear model, though
this has not been rigorously demonstrated.
Now, let us consider the non-linear analogue of the generalized least squares
estimator, which I will call the non-linear generalized least squares (NLGLS)

on an approximation of A by a circular matrix. [See Amemiya and Fuller (1967,
p. 527).]
Hannan proves the strong consistency of his non-linear spectral estimator
obtained by minimizing the right-hand side of (3.7) under the assumptions (2.6),
356
T. Amemiya
(2.12), and the new assumption
f
CfAcAf(r+s)(c*)
converges uniformly in
c, ,
c2 E B
for every integer S.
f
(3.8)
Note that this is a generalization of the assumption (2.11). However, the assump-
tion (3.8) is merely sufficient and not necessary. Hannan shows that in the model
y, = OL, + (Y~COS& + cr,sin&t + u,,
(3.9)
assumption (3.8) does not hold and yet b is strongly consistent if we assume (3.4)
and 0 < /?a < T. In fact, T(fi - &) converges to zero almost surely in this case.
In proving the asymptotic normality of his estimator, Hannan needs to gener-
alize (2.20) and (2.21) as follows:
+c
$i,,
*I,,
converges uniformly in c, and cZ
in an open neighborhood of &
(3.10)
and

(3.14)
we can write
A
= (2?r)-‘/Y,g(w)+(o)*dF(w) and
B =
(2a)-‘/l,+(o)dF(w).
Ch. 6: Non -linear Regression Models
357
In the model (3.9), assumptions (3.10) and (3.11) are not satisfied; nevertheless,
Hannan shows that the asymptotic normality holds if one assumes (3.4) and
0 < & < 7r. In fact, J?;T(b - &) , normal in this case.
An interesting practical case is where I#B(W) = a)‘, where g(w) is a con-
sistent estimator of g(o). I will denote this estimator by b(e). Harman proves
that B(2) and b(Z) have the same asymptotic distribution if g(w) is a rational
spectral density.
Gallant and Goebel (1976) propose a NLGLS estimator of the autocorrelated
model which is constructed in the time domain, unlike Hannan’s spectral estima-
tor. In their method, they try to take account of the autocorrelation of {u,} by
fitting the least squares residuals ti, to an autoregressive model of a finite order.
Thus, their estimator is a non-linear analogue of the generalized least squares
estimator analyzed in Amemiya (1973a).
The Gallant-Goebel estimator is calculated in the following steps. (1) Obtain
the NLLS estimator 8. (2) Calculate li = y - f(b). (3) Assume that (u,} follow an
autoregressive model of a finite order and estimate the coefficients by the least
squares regression of z?, on
zi,_ ,
, zi,_ 2,. . . .
(4) Let 2 be the variance-covariance
matrix of u obtained under the assumption of an autoregressive model. Then we
can find a lower triangular matrix