HOMOTOPY CONTINUATION METHODS:
An Algorithm for the Fixed Point and Newton Homotopy Methods with Some Examples
Francesco Decarolis Ricardo Mayer Martin Santamaria
fdc@uchicago.edu ricardom@uchicago.edu msantama@uchicago.edu
ABSTRACT
Homotopy Continuation Methods provide a useful
approach to find the zeros of f : R
n
→R
n
in a globally
convergent way. Homotopy Methods transform a hard
problem into a simpler one whit easily calculated zeros
and then gradually deform this simpler problem into
the original one computing the zeros of the intervening
problems and eventually ending with a zero of the
original problem.
This paper presents a description of a Matlab code that
can implement the Fixed Point Homotopy and the
Newton Homotopy methods to find the zeros of any
function f : R
n
→R
n
.
Four examples are also provided for illustrative
purposes. Two univariate examples, one of which
regarding a General Equilibrium problem, and two
multivariate ones show how the Matlab code can be
easily applied to different problems.
You will need Matlab's Symbolic Math Toolbox to run
this code.
1. INTRODUCTION: the Homotopy Continuation Methods in General
This paper accompanies a code that can implement the Fixed Point Homotopy and the
Newton Homotopy methods to find the zeros of any function f : R
n
→R
n
. Homotopy
methods belong to the family of continuation methods and like all these methods they
represent a way to find a solution to a problem by constructing a new problem, simpler
than the original one, and then gradually deforming this simpler problem into the original
one keeping track of the series of zeros that connect the solution of the simpler problem
to that of the original, harder one. The greatest advantage of the Homotopy Methods is
that, under some conditions, they offer a way to have a globally convergent method to
find the zeros of any function f : R
n
→R
n
.
Judd (1998), chapter 5, contains a detailed description of the computational issues
involved in the use of these Homotopy Continuation Methods. In this paper we follow the
discussion in Judd (1998) to construct a simple code that allows to use the Fixed Point
Homotopy (FPH) and the Newton Homotopy (NH) to find the zeros of f : R
n
→R
n
.
All the Homotopy methods are based on the construction of a function, H(x,t),
H: R
n+1
→R
n
, H∈C
0
(R
n+1
). This function H , which is called the Homotopy function,
serves the objective of continuously deforming a function g into f and can be any
continuous function such that: H(x,0)=g(x) and H(x,1)=f(x). Hence we construct H(x,0)
in such a way that its zeros are easily found while we also require that, once the
parameter t is equal to 1, then H coincides with the original function f .
Among the various Homotopy functions that are generally used the Fixed Point
Homotopy and the Newton Homotopy are simple and powerful methods that can be
successfully applied to several different problems. In the next sections we describe how
these two Homotopy functions are constructed and we illustrate how the Matlab code that
we wrote applies these methods in an algorithm to find the zeros of f : R
n
→R
n
. Finally
we illustrate how to actually use the Matlab code we wrote through four examples, both
univariate and multivariate.
2. FIXED POINT HOMOTOPY, NEWTON HOMOTOPY AND THEIR
IMPLEMENTATION
The code described in the next section is an algorithm that finds zeros of f : R
n
→R
n
through an Homotopy method. The two Homotopy methods between which our code
allows to chose are the following:
(a) Fixed Point Homotopy: H(x,t)=(1-t) (x-x
0
)+tf(x) for some x
0
.
(b) Newton Homotopy: H(x,t)= f(x) - (1-t) f(x
0
) for some x
0
.
The first Homotopy function gradually deforms the function (x-x
0
) into f(x), while the
second one is a function that for t=0 has its zero at x=x
0
while at t=1 it is f(x). The idea of
our code will be to follow the path connecting the zeros of H(x,t) from t=0, where they
are as described above, until t=1, where we will have found the zeros of the original
problem. If H(x,0) and H(x,1) have zeros we will try to develop a way to follow a
continuous path, which we indicate as H
-1
(0), H
-1
(0)={(x,t)│H(x,t)=0}, that connects one
to the other.
The procedure we follow in order to create a code that traces the path H
-1
(0) is based on
the parametric approach described in Judd (1998). Hence we view x and t as both
functions of the same parameter s and we require that a parametric path satisfies
H(x(s),t(s))=0 for all s. Being on the path thus implies that we have an y(s)=(x(s),t(s))
such that H is equal to zero hence for changes in x and t it must be that y obeys the
following system differential equations (see Garcia and Zangwill (1981)):
dy
i
ds
−1
i
det
∂H
∂y
y
−i
, i 1,...,n1
where ( · )
-i
means that we remove the ith column.
For the system defined by the above equation to arrive to t(s)=1 it must be that we never
find any singular matrix (((∂H)/(∂y))(y)
-i
) before that. The condition that insures that this
is the case is that the Homotopy function H is regular, where we define H to be regular
whenever H
y
=(H
x
,H
t
) is of full rank (i.e. rank n) at all points in H
-1
(0) and y≡(x,t).
The code that we present in the next section builds on this parametric approach to
describe a path in the (x,t) space that connects the zeros of H(x,0) to those of H(x,1).
3. THE CODE
The code involves using only 2 m-files, only one of them written by the user. These are:
• A User's Script (US) : this is the only file a user will need to run. Here the user must
write the equations he wants to solve (it may be only one equation) as long with the
other inputs needed by htopy.m. This program calls the function htopy, get's the
results and further calls fsolve to refine the solution delivered by htopy. We suggest
the user first read carefully and understand well the example file called
'examplebook.m' . This example show how to solve an univariate problem and
contains extensive comments. Later, the reader should take a look at any other
example just to see the small changes needed to solve a multivariate problem. Note
that there is the freedom to choose any letter (like x, q, z, etc., or even multi-characters
expressions like x1, x2, xx1, xx2, etc.) for the independent variables
with the
exception of y and t , because those will be names used internally by htopy. The order
in which the user inputs these variables will be the one in which the plots and of the
elements in the Homotopy solution vector, hopt, will appear.
• The m-file function htopy.m. The function is relatively well-documented and we
encourage the user to have a glance of at least the first comments included in the file:
it provides a description of the inputs needed to call the function, the output it gives,
describes the organization of the whole algorithm and the choices made during the
execution of the program. Is short, htopy has two main parts:
(a) The first uses the Symbolic Math Toolbox's capabilities to construct the system of
ODEs defining the solution's path. This is, the user only have to know his initial
system of equation and not the possibly quite different system of ODEs: this is
computed on-the-fly by htopy. In particuar it will compute a symbolic Jacobian of the
system and symbolic Determinants of parts of the Jacobian. It never writes a new m-
file: it will store the system of ODEs as an inline function to feed the ODE-solver
later. Some work inside the function is spent in converting symbolic object into string
objects suitable for inline functions.
(b) The second part takes such system of ODEs -now in form an inline function- and
calls the ODE-solver specified by the user. It does the path-tracking and explores
whether we have reached a solution and whether the solution is an acceptable one.
Conditional on success, returns the value of the candidate solution. If something fails
it gives an appropriate warning and information to the User so she can decide how to
proceed next.
4. FOUR APPLICATIONS OF THE HOMOTOPY CODE
Here we present four examples that will illustrate how to use the code described above.
The first two examples are univariate. In particular the first one is the same used in Judd
(1998), chapter 5, to present the Homotopy methods. Here it will be particularly useful to
exemplify the differences between the Fixed Point and the Newton Homotopy. The
second one, instead, is a simple General Equilibrium and we show how the code is used
to find the equilibrium prices. Then, finally, the third and fourth examples illustrate how
to apply the htopy.m code to multivariate problems.
4.1 Example1 (Taken from Judd (1998)).
The problem is to find the value of x for which f(x)=0, where f(x) is given by:
f(x)=2*x - 4 + sin(2*pi*x)
The unique solution of this problem is x=2. Below we show the results that we get from
the code reported in the Appendix at the end of this paper. We pick as our initial guess
for x, x0=0. We further choose to have 70 checks and UpS equal to 2.5. For the two
Homotopy methods that our htopy.m function can compute we get the following results:
(a) Fixed Point Homotopy
Computed Solutions
Homotopy Solution: 2.0008e+000
Gauss-Newton Solution: 2.0000e+000
评论5