Notes on the fixed point framework

2018-03-27

In this vignette you can find some notes and examples on the use of the fixed point function and associated tools. The function itself is simple enough; it evaluates:

\[ x_{n + 1} = f(x_{n}) \]

until a convergence criterion is reached. A convergence criterion is defined as a function:

\[ convCrit(x_{n + 1}, x_{n}) = TRUE \]

if convergence is reached and FALSE otherwise. As a simple example consider the function

\[ f(x) = 1 + \frac{1}{x} \]

To find the value of \(x\) for which \(f(x) = x\) we can use fixedPoint.

library("saeRobust")
## Loading required package: aoos
library("magrittr")
convCrit <- function(xn1, xn0) abs(xn0 - xn1) < 0.001
fixedPoint(function(x) 1 + 1 / x, rnorm(1), convCrit = convCrit)
## [1] 1.618205

Square Root

Another example which can be used is to find the square root. The square root of \(p\) is the fixed point of the function

\[ f(x) = \frac{p}{x} \]

In this example, however, the function does not necessarily converge. To show this we modify our convCrit function to stop the algorithm after 10 iterations:

sqrtFp <- function(p) {
  force(p)
  function(x) p / x
}
fixedPoint(sqrtFp(2), 2, addMaxIter(convCrit, 10))
## [1] 2

Surprisingly we are not even close to the actual value. To check the values in each iteration you can modify your fp function such that it saves the history of values in all iterations:

fixedPoint(addHistory(sqrtFp(2)), 2, addMaxIter(convCrit, 10))
## [1] 2
## attr(,"history")
##     [,1]
## res    1
## res    2
## res    1
## res    2
## res    1
## res    2
## res    1
## res    2
## res    1
## res    2

What happens is that the functions oscillates between 1 and 2. To overcome this average damping can be applied, which can also be beneficial with regards to the speed of convergence:

fixedPoint(addHistory(addAverageDamp(sqrtFp(2))), 2, addMaxIter(convCrit, 10))
## [1] 1.414214
## attr(,"history")
##         [,1]
## res 1.500000
## res 1.416667
## res 1.414216
## res 1.414214

Square root solved with Newton-Raphson

The Newton Raphson is frequently used to fit statistical models. It is a special case of the fixed point algorithm and exists as a convenience in this package. The algorithm is defined as:

\[ x_{n + 1} = x_n + (f'(x_n))^{-1} f(x_n) \]

To find the square root it is sufficient to find the root of

\[ f(x) = x^2 - p \]

So we can implement the algorithm as follows:

sqrtNr <- function(.p) {
  force(.p)
  f <- function(x) x^2 - .p
  f1 <- function(x) 2 * x
  as.list(environment())
}

newtonRaphson(sqrtNr(2), 2, convCrit = convCrit)
## [1] 1.414214