Skip to content

Commit

Permalink
updating documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
@s-baumann
s-baumann committed Jan 8, 2019
1 parent c8b76b8 commit 92efa97
Show file tree
Hide file tree
Showing 5 changed files with 69 additions and 103 deletions.
2 changes: 2 additions & 0 deletions .travis.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ julia:
notifications:
email: false

# This does the documentation
jobs:
include:
- stage: "Documentation"
Expand All @@ -22,3 +23,4 @@ jobs:
- julia --project=docs/ docs/make.jl

coveralls: true
coverage: true
1 change: 1 addition & 0 deletions docs/make.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ makedocs(
sitename = "FixedPointAcceleration",
modules = [FixedPointAcceleration],
pages = ["index.md",
"1_FixedPoints.md",
"2_Algorithms.md",
"3_UsingAdvice.md",
"99_refs.md"]
Expand Down
10 changes: 10 additions & 0 deletions docs/src/1_FixedPoints.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@

# 1 Fixed point acceleration

A fixed point problem is one where we look for a vector, $\hat{X} \in \Re^N$, so that for a given function $f: \Re^N \rightarrow \Re^N$ we have:

$$f(\hat{X}) = \hat{X}$$

If $f: \Re^1 \rightarrow \Re^1$ and thus any solution $\hat{X}$ will be a scalar then one way to solve this problem would be to use a rootfinder on the function $g(x) = f(x) - x$ or to use an optimiser to minimise $h(x) = (f(x) - x)^2$. These techniques will not generally work however if $f : N^a \rightarrow N^a$ where $a$ is large. Consider for instance using a multidimensional Newtonian optimiser to minimise $h(x) = (f(x) - x)^2$. The estimation of gradients for each individual dimension may take an infeasibly long time. In addition this method may not make use all available information. Consider for instance that we know that the solution for $x$ will be an increasing vector (so $x_i > x_j$ for any entries of $x$ with $i > j$) but has many entries. This information can be preserved and used in the vector acceleration algorithms that we present but would be more difficult to exploit in a standard optimisation algorithm.

**FixedPointAcceleration.jl** implements eight algorithms for finding fixed points. The first algorithm implemented in this package is the "simple" method which merely takes the output of a function and feeds it back into the function. For instance starting with a guess of $x_0$, the next guess will be $x_1 = f(x_0)$. The guess after that will be $x_2 = f(x_1)$ and so on. In some conditions $f$ will be a contraction mapping and so the simple method will be guaranteed to converge to a unique fixed point (Stokey, Lucas & Prescott 1989). Even when this is the case however the simple method may only converge slowly which can be inconvenient. The other seven methods this package implements are designed to be faster than the simple method but may not be convergent for every problem.
149 changes: 55 additions & 94 deletions docs/src/3_UsingAdvice.md
View file
Original file line number Diff line number Diff line change
@@ -1,150 +1,111 @@

# 3. Using the FixedPoint package
# 3. Using the FixedPointAcceleration package

## 3.1 Basic examples of using FixedPoint
## 3.1 Basic examples of using FixedPointAcceleration

### The Babylonian method for finding square roots.

Now we will demonstrate how **FixedPoint** can be used for simple problems. For the simplest possible case consider we want to estimate a square root using the Babylonian method. To find the square root of a number $x$, given an initial guess $t_0$ the following sequence converges to the square root:
Now we will demonstrate how **FixedPointAcceleration** can be used for simple problems. For the simplest possible case consider we want to estimate a square root using the Babylonian method. To find the square root of a number $x$, given an initial guess $t_0$ the following sequence converges to the square root:

$t_{n+1} = \frac{1}{2} \left[ t_n + \frac{x}{t_n} \right]$

This is a fast converging and inexpensive sequence which probably makes an acceleration algorithm overkill but for sake of exposition we can implement this in **FixedPoint**. In the next code block we find the square root of 100 with the simple method:
This is a fast converging and inexpensive sequence which probably makes an acceleration algorithm overkill but for sake of exposition we can implement this in **FixedPointAcceleration**. In the next code block we find the square root of 100 with the simple method:

```
using FixedPointAcceleration
SequenceFunction = function(tn){0.5*(tn + 100/tn)}
SequenceFunction(x) = 0.5 .* (x .+ 100 ./ x)
Initial_Guess = 6.0
FP_Simple = fixed_point(SequenceFunction,Initial_Guess; Method = "Simple")
FP_Simple = fixed_point(SequenceFunction, Initial_Guess; Algorithm = Simple)
```

We can also solve for a vector of fixed points at the same time. For instance every square root from 1 to 100.

```r
NumbersVector = 1:100
SequenceFunction = function(tn){0.5*(tn + NumbersVector/tn)}
FP_SEA = FixedPoint(Function = SequenceFunction, Inputs = 1:100, Method = "RRE")
```
NumbersVector = collect(1:100)
SequenceFunction(x) = 0.5 .* (x .+ NumbersVector ./ x)
Initial_Guess = repeat([10],100)
FP_SEA = fixed_point(SequenceFunction, Initial_Guess; Algorithm = RRE)
```
Note that in this case the RRE method is being applied elementwise.

### Vectorised functions

The utility of the acceleration algorithms contained in **FixedPoint** are more apparent when applied to vectorised functions with cross dependency. For a simple example consider the below function where each entry of the vector depends on both entries of the previous iterate.

```r
SimpleVectorFunction = function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])}
FP_Simple = FixedPoint(Function = SimpleVectorFunction, Inputs = c(0.3,900), Method = "Simple")
FP_Anderson = FixedPoint(Function = SimpleVectorFunction, Inputs = c(0.3,900), Method = "Anderson")
```
SimpleVectorFunction(x) = [0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2]]
Initial_Guess = [0.3,900]
FP_Simple = fixed_point(SimpleVectorFunction , Initial_Guess; Algorithm = Simple)
FP_Anderson = fixed_point(SimpleVectorFunction, Initial_Guess; Algorithm = Anderson)
```
This function takes 105 iterates to find a fixed point with the simple method but only 14 with the Anderson acceleration method.

## 3.2 Easily changing algorithm

A key priority in writing **FixedPoint** is for ease in changing fixed point method. The algorithm can be stopped and then restarted without having to repeat earlier iterations. This can be seen in the simple example below:

```r
Run1 = FixedPoint(Function = function(x){x^0.2}, Inputs = c(0.3, 0.4,0.5),
Method = "Aitken")
cat(paste("Running this all at once takes ", length(Run1$Convergence), " iterations and
convergences to ", Run1$Convergence[length(Run1$Convergence)], "\n"))
We can "chain" together different calls to the fixed_point function in order to switch acceleration algorithm
at any point. For instance consider the following function and initial guess at a fixed point:
```

func(x) = [0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2]]
Initial_Guess = [1.1,2.2]
```
## Running this all at once takes 9 iterations and
## convergences to 1.11022302462516e-15
Now we can initially do two simple iterates. Then do three iterates with the MPE method. Then one with the simple method and then finish with the RRE method. This can be done in the following way:
```

```r
Run2_1 = FixedPoint(Function = function(x){x^0.2}, Inputs = c(0.3, 0.4,0.5),
Method = "Aitken", MaxIter = 4)
# Check some condition about the convergence or change algorithm.
Run2_2 = FixedPoint(Function = function(x){x^0.2}, Inputs = Run2_1$Inputs,
Outputs = Run2_1$Outputs, Method = "Aitken")
cat(paste("Running this with a break takes ", length(Run2_2$Convergence), " iterations and
convergences to ", Run2_2$Convergence[length(Run2_2$Convergence)], "\n"))
fp_chain = fixed_point(func, Initial_Guess; Algorithm = Simple, MaxIter = 2)
fp_chain = fixed_point(func, fp_chain; Algorithm = MPE, MaxIter = 3)
fp_chain = fixed_point(func, fp_chain; Algorithm = Simple, MaxIter = 1)
fp_chain = fixed_point(func, fp_chain; Algorithm = RRE, MaxIter = 100)
```

Now as it turns out The MPE (and RRE) does simple iterates except for every iterate that is a multiple of the ExtrapolationPeriod (7 by default). And so there is no difference from the above sequence of iterates and just doing all iterates with the RRE. This can be verified with the following:
```
## Running this with a break takes 9 iterations and
## convergences to 1.11022302462516e-15
fp_nochain = fixed_point(func, Inputs; Algorithm = RRE, MaxIter = 100)
fp_chain.Iterations_ == fp_nochain.Iterations_
all(abs.(fp_chain.Inputs_ .- fp_nochain.Inputs_) .< 1e-14)
```
This does highlight that there is no waste in changing fixed_point algorithm in this way. No iterates are reevaluated.

This can be useful in problems where it is desired to switch extrapolation methods half way through. For instance in the consumption smoothing case presented later, the Vector Epsilon Algorithm converges quickly to a convergence level around $10^{-4}$ but then spends around 100 iterations hovering around this level without reducing further. In this case the algorithm could be changed to the simple method or Anderson method once the algorithm hit some target level of convergence.

Another example is parallelisation. The Anderson method in particular is well suited to use the results of previous iterates and these iterates do not need to be generated sequentially. If a suitable collection of guess vectors can be generated[^4] then all could be evaluated in parallel and recombined before again being used by the Anderson algorithm to generate a new guess. An example of this is presented later when the package is used to find the value function for a consumer in a consumption smoothing problem.
Changing algorithms can be useful in some cases where an error occurs. For instance consider we are trying to find the
fixed point of the following function:
```
simple_vector_function(x) = [0.5*sqrt(x[1] + x[2]), 1.5*x[1] + 0.5*x[2]]
Inputs = [0.3,900]
fp = fixed_point(simple_vector_function, Inputs; Algorithm = Anderson)
```
Inspecting this fp object reveals an error after the 3rditeration because Anderson tries to use a negative value for both x entries which results in the square root of a negative number. We can switch to simple iterations to get closer to the fixed point at which point Anderson will no longer try negative numbers. This will fix this.
```
fp = fixed_point(simple_vector_function, fp; Algorithm = Simple, MaxIter = 7)
fp = fixed_point(simple_vector_function, fp; Algorithm = Anderson)
```

[^4]: For instance a set of guesses could be done by using different dampening values or by running the Anderson algorithm several times using different subsets of previous iterates to generate new guesses.
## 3.3 Graceful error handling

## 3.4 Graceful error handling
Hopefully **FixedPointAcceleration** is well tested enough that most kind of errors will be rare. **FixedPointAcceleration** also offers an option (ReplaceInvalids) to ensure that no acceleration algorithm generates guess vectors that contain NAs or Infs. This option can be set to ReplaceVector which will replace an extrapolated vector containing missings, NANs or Infs by the vector output in the previous iterate. If it is set to ReplaceElement then it will replace the individual elements that are missings, NANs or Infs by the corresponding elements in the output of the previous iterate.

Hopefully **FixedPoint** is well tested enough that most kind of errors will be rare. **FixedPoint** also offers an option (``ReplaceInvalids'') to ensure that no acceleration algorithm generates guess vectors that contain NAs or Infs. This option can be set to ``ReplaceVector'' which will replace an extrapolated vector containing NAs or Infs by the vector output in the previous iterate. If it is set to ``ReplaceElement'' then it will replace the individual elements that are NA or Inf by the corresponding elements in the output of the previous iterate.[^15]
Errors are likely however in cases where inputs functions have a restricted domain. For example this may include functions that require the input vector to have a particular shape (ie concavity) or functions where the input vector must be strictly positive. For a simple example consider the vectorised function we introduced in section 3.1. Now rather than

[^15]: With these options being set the **FixedPoint** should never throw an error when applied to a function which is guaranteed to map: $\Re^N \rightarrow \Re^N$ for some $N$ unless the error comes from underflow or overflow. If you do encounter such a testcase it likely reflects a bug, we would much appreciate it if you could send a minimal working example to the package maintainer.
$x^\prime[1] = \frac{\sqrt{\vert x[1] + x[2] \vert}}{2}$ we have $x^\prime[1] = \frac{\sqrt{ x[1] + x[2] }}{2}$

Errors are likely however in cases where inputs functions have a restricted domain. For example this may include functions that require the input vector to have a particular shape (ie concavity) or functions where the input vector must be strictly positive. For a simple example consider the vectorised function we introduced in section 3.1. Now rather than $x^\prime[1] = \frac{\sqrt{\vert x[1] + x[2] \vert}}{2}$ we have $x^\prime[1] = \frac{\sqrt{ x[1] + x[2] }}{2}$ where the output x has a prime and the inputs has no prime symbol. $x^\prime[1]$ here is no longer real valued if the sum of the previous iterate is negative. This is what occurs in the 5th iterate of the Anderson method applied to this problem.
where the output x has a prime and the inputs has no prime symbol. $x^\prime[1]$ here is no longer real valued if the sum of the previous iterate is negative. This is what occurs in the 5th iterate of the Anderson method applied to this problem.

The FixedPoint function handles these situations gracefully by saving all previous results as well as the proposed new vector that lead to the error. In these cases the "Finish" message in the output describes the form of error, the "NewInputVector" list entry contains the attempted vector that was input to the function and the "NewOutputVector" list entry contains the result of applying the function to NewInputVector. This output information can be seen in the example below.
The FixedPoint function handles these situations gracefully by saving all previous results as well as the proposed new vector that lead to the error. In the event of such an error the FailedEvaluation_ member of the returned FixedPointResults struct will describe the issue.

This information is useful in order to diagnose the issue. In this case we might decide to modify the function to insert the absolute value function with the reasoning that the same fixed point will likely result (which we could later verify). This also allows a user to run one method until an error occurs and then switch methods. This is demonstrated below.


```r
SimpleVectorFunction = function(x){c(0.5*sqrt(x[1] + x[2]), 1.5*x[1] + 0.5*x[2])}
FPSolution = FixedPoint(Function = SimpleVectorFunction, Inputs = c(0.3,900),
Method = "Anderson")
```

```
## Warning in sqrt(x[1] + x[2]): NaNs produced
```

```r
# We can see the output of the FixedPoint function in cases like this where it ends due
# to an error
FPSolution
SimpleVectorFunction(x) = [0.5*sqrt(x[1] + x[2]), 1.5*x[1] + 0.5*x[2]]
Initial_Guess = [0.3,900]
FPSolution = FixedPoint(SimpleVectorFunction, Initial_Guess; Algorithm = Anderson)
```

```
## $Inputs
## [,1] [,2] [,3] [,4]
## [1,] 0.3 15.0025 7.350817 3.191655
## [2,] 900.0 450.4500 82.469248 9.574966
##
## $Outputs
## [,1] [,2] [,3] [,4]
## [1,] 15.0025 10.78717 4.738672 1.786520
## [2,] 450.4500 247.72875 52.260849 9.574966
##
## $Convergence
## [1] 449.550000 202.721250 30.208399 1.405135
##
## $FixedPoint
## [1] NA
##
## $Finish
## [1] "New output vector contains NAs"
##
## $NewInputVector
## [1] -1.085113 -3.255338
##
## $NewOutputVector
## [1] NaN -3.255338
```

```r
# We can use this information to decide to switch to the simple method.
# No error results as the simple method doesn't extrapolate.
FPSolution = FixedPoint(Function = SimpleVectorFunction, Inputs = FPSolution$Inputs,
Outputs = FPSolution$Outputs, Method = "Simple", MaxIter = 5)
FPSolution = FixedPoint(SimpleVectorFunction, FPSolution; Algorithm = Simple, MaxIter = 5)
# Now we switch to the Anderson Method again. No error results because we are
# close to fixed point.
FPSolution = FixedPoint(Function = SimpleVectorFunction, Inputs = FPSolution$Inputs,
Outputs = FPSolution$Outputs, Method = "Anderson")
FPSolution = FixedPoint(SimpleVectorFunction, FPSolution; Algorithm = Anderson)
```



## 3.5 Convergence by constant increments
## 3.4 Convergence by constant increments

Most of the methods included in this function will fail in finding the fixed point of a function that converges by a fixed increment.
For instance we may have a function that takes $x$ and returns $x$ shifted 1 unit (in Euclidian norm) in a straight line
Expand All @@ -165,4 +126,4 @@ In these cases it may be possible to change the function to make it converge by
same set of fixed points. This is shown in the perceptron example in section 4.3. In other cases where it is not possible
to modify the function, it is advisable to use the simple method.

[^5]: When these complications arise the ReplaceInvalids method can be used to revert to a simple iterate or to change individual elements to the corresponding values in a simple iterate. This is as described in section 3.4.
[^5]: When these complications arise the ReplaceInvalids method can be used to revert to a simple iterate or to change individual elements to the corresponding values in a simple iterate. This is as described in section 3.3.
10 changes: 1 addition & 9 deletions docs/src/index.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,18 +8,10 @@ In this paper section 1 starts by with a brief explanation of fixed points befor

Section 4 then presents several applications of these fixed point algorithms in economics, asset pricing and machine learning. Finally section 5 presents a convergence comparison showing how many iterations each algorithm takes in bringing each problem to its fixed point for each of the applications presented in section 4.

# 1 Fixed point acceleration

A fixed point problem is one where we look for a vector, $\hat{X} \in \Re^N$, so that for a given function $f: \Re^N \rightarrow \Re^N$ we have:

$$f(\hat{X}) = \hat{X}$$

If $f: \Re^1 \rightarrow \Re^1$ and thus any solution $\hat{X}$ will be a scalar then one way to solve this problem would be to use a rootfinder on the function $g(x) = f(x) - x$ or to use an optimiser to minimise $h(x) = (f(x) - x)^2$. These techniques will not generally work however if $f : N^a \rightarrow N^a$ where $a$ is large. Consider for instance using a multidimensional Newtonian optimiser to minimise $h(x) = (f(x) - x)^2$. The estimation of gradients for each individual dimension may take an infeasibly long time. In addition this method may not make use all available information. Consider for instance that we know that the solution for $x$ will be an increasing vector (so $x_i > x_j$ for any entries of $x$ with $i > j$) but has many entries. This information can be preserved and used in the vector acceleration algorithms that we present but would be more difficult to exploit in a standard optimisation algorithm.

**FixedPointAcceleration.jl** implements eight algorithms for finding fixed points. The first algorithm implemented in this package is the "simple" method which merely takes the output of a function and feeds it back into the function. For instance starting with a guess of $x_0$, the next guess will be $x_1 = f(x_0)$. The guess after that will be $x_2 = f(x_1)$ and so on. In some conditions $f$ will be a contraction mapping and so the simple method will be guaranteed to converge to a unique fixed point (Stokey, Lucas & Prescott 1989). Even when this is the case however the simple method may only converge slowly which can be inconvenient. The other seven methods this package implements are designed to be faster than the simple method but may not be convergent for every problem.

```@contents
pages = ["index.md",
"1_FixedPoints.md",
"2_Algorithms.md",
"3_UsingAdvice.md",
"99_refs.md"]
Expand Down

0 comments on commit 92efa97

Please sign in to comment.