The fixed point problem is to find a solution to the equations \( x=f(x) \), where \( f \) is a continuous function, \( f:\mathbb{R}^n\Rightarrow\mathbb{R}^n \). Many systems of equations \( g(x)=0 \) can be rewritten as \( x-f(x)=0 \), i.e. as \( x=f(x) \).
We can try to solve the fixed problem iteratively by using iterations of the form \( x^{(k+1)}=f(x^{(k)}). \) If the sequence \( x^{(k)} \) converges to, say, \( x_\infty \), then \( x_\infty=f(x_\infty) \), i.e. \( x_\infty \) is a fixed point of \( f \), a solution of the fixed point problem.
Below is the code for a general fixed point iterator. As arguments it needs an initial estimate xinit and the function f being iterated. The function f may require additional arguments. In addition some parameters, with default values, are used to determine the precision, bound the number of iterationss, and regulate the amount of output.
myIterator <- function(xinit, f, eps = 1e-06, itmax = 100, verbose = FALSE,
...) {
xold <- xinit
itel <- 1
repeat {
xnew <- f(xold, ...)
if (verbose)
cat("Iteration: ", formatC(itel, width = 3, format = "d"), "xold: ",
formatC(xold, digits = 8, width = 12, format = "f"), "xnew: ",
formatC(xnew, digits = 8, width = 12, format = "f"), "\n")
if ((supDist(xold, xnew) < eps) || (itel == itmax)) {
return(xnew)
}
xold <- xnew
itel <- itel + 1
}
}
The function uses a small subroutine to measure distance between successive iterations, which can easily be replaced by other ways of measuring distance, and thus with another way of establishing convergence.
supDist <- function(x, y) return(max(abs(x - y)))
We can use fixed point iterations to rapifdly compute the square root of any number. The function we use for the iterations is
fSqrt <- function(x, y) (x + (y/x))/2
Let's use this to compute \( \sqrt{12} \).
sy <- myIterator(2, fSqrt, eps = 1e-06, itmax = 100, verbose = TRUE, y = 12)
## Iteration: 1 xold: 2.00000000 xnew: 4.00000000
## Iteration: 2 xold: 4.00000000 xnew: 3.50000000
## Iteration: 3 xold: 3.50000000 xnew: 3.46428571
## Iteration: 4 xold: 3.46428571 xnew: 3.46410162
## Iteration: 5 xold: 3.46410162 xnew: 3.46410162
In Multidimensional Scaling we minimize a function of the form \[ \sigma(x)=\frac12\sum_{i=1}^nw_i(\delta_i-\sqrt{x'A_ix})^2, \] where the \( \delta_i \) and \( w_i \) are known positive numbers and the \( A_i \) are known positive semi-definite matrices which satisfy \( \sum_{i=1}^n w_iA_i=I \).
Define \( d_i(x)=\sqrt{x'A_ix} \) and \[ B(x)=\sum_{i=1}^n w_i\frac{\delta_i}{d_i(x)}A_i. \] By differentiating the loss function we see that local minima must satisfy the stationary equations \( x=B(x)x \), which shows that we are looking for a fixed point of \( f(x)=B(x)x \). The code for the function is
mds <- function(x, a, delta) {
d <- apply(a, 3, function(z) sqrt(sum(x * z %*% x)))
b <- sliceSums(outer(a, delta/d)[, , 1, ])
return(drop(b %*% x))
}
This uses the subroutine sliceSums(), similar to rowSums() and columnSums().
sliceSums <- function(a) {
n <- dim(a)
b <- matrix(0, n[1], n[1])
for (i in 1:n[3]) {
b <- b + a[, , i]
}
return(b)
}
Make an array of three \( 4\times 4 \) matrices \( A_i \) that add up to \( I \).
set.seed(12345)
a <- array(0, c(4, 4, 3))
a[, , 1] <- crossprod(matrix(rnorm(16), 4, 4))
a[, , 2] <- crossprod(matrix(rnorm(16), 4, 4))
a[, , 3] <- crossprod(matrix(rnorm(16), 4, 4))
aa <- sliceSums(a)
hh <- eigen(aa)
hv <- hh$vectors
hw <- hh$values
hz <- sqrt(outer(hw, hw))
a[, , 1] <- (t(hv) %*% a[, , 1] %*% hv)/hz
a[, , 2] <- (t(hv) %*% a[, , 2] %*% hv)/hz
a[, , 3] <- (t(hv) %*% a[, , 3] %*% hv)/hz
sliceSums(a)
## [,1] [,2] [,3] [,4]
## [1,] 1.000e+00 1.527e-16 -4.580e-16 1.388e-17
## [2,] 1.579e-16 1.000e+00 -1.735e-16 3.539e-16
## [3,] -3.886e-16 -1.180e-16 1.000e+00 -2.359e-16
## [4,] 9.714e-17 2.984e-16 -2.220e-16 1.000e+00
Now run our fixed point iterations.
myIterator(1:4, mds, eps = 0.001, verbose = TRUE, a = a, delta = c(1, 2, 3))
## Iteration: 1 xold: 1.00000000 2.00000000 3.00000000 4.00000000 xnew: 1.35773413 2.38493320 0.40548519 2.41719162
## Iteration: 2 xold: 1.35773413 2.38493320 0.40548519 2.41719162 xnew: 1.81471484 2.46733696 0.12552570 3.00393225
## Iteration: 3 xold: 1.81471484 2.46733696 0.12552570 3.00393225 xnew: 2.03235214 2.28196303 -0.14645695 3.31689866
## Iteration: 4 xold: 2.03235214 2.28196303 -0.14645695 3.31689866 xnew: 2.18193053 2.10100837 -0.37755080 3.55116991
## Iteration: 5 xold: 2.18193053 2.10100837 -0.37755080 3.55116991 xnew: 2.27319040 1.93173883 -0.55906455 3.69871107
## Iteration: 6 xold: 2.27319040 1.93173883 -0.55906455 3.69871107 xnew: 2.32447733 1.78630465 -0.69432926 3.78299071
## Iteration: 7 xold: 2.32447733 1.78630465 -0.69432926 3.78299071 xnew: 2.35166773 1.66884163 -0.79243393 3.82822962
## Iteration: 8 xold: 2.35166773 1.66884163 -0.79243393 3.82822962 xnew: 2.36546458 1.57761148 -0.86283245 3.85150227
## Iteration: 9 xold: 2.36546458 1.57761148 -0.86283245 3.85150227 xnew: 2.37217525 1.50840643 -0.91324715 3.86304393
## Iteration: 10 xold: 2.37217525 1.50840643 -0.91324715 3.86304393 xnew: 2.37526061 1.45664155 -0.94941456 3.86852231
## Iteration: 11 xold: 2.37526061 1.45664155 -0.94941456 3.86852231 xnew: 2.37654659 1.41824772 -0.97543957 3.87094909
## Iteration: 12 xold: 2.37654659 1.41824772 -0.97543957 3.87094909 xnew: 2.37697272 1.38991844 -0.99422490 3.87188593
## Iteration: 13 xold: 2.37697272 1.38991844 -0.99422490 3.87188593 xnew: 2.37701205 1.36908335 -1.00782184 3.87212722
## Iteration: 14 xold: 2.37701205 1.36908335 -1.00782184 3.87212722 xnew: 2.37689660 1.35379203 -1.01768563 3.87206852
## Iteration: 15 xold: 2.37689660 1.35379203 -1.01768563 3.87206852 xnew: 2.37673572 1.34258488 -1.02485396 3.87190091
## Iteration: 16 xold: 2.37673572 1.34258488 -1.02485396 3.87190091 xnew: 2.37657753 1.33437869 -1.03007054 3.87171217
## Iteration: 17 xold: 2.37657753 1.33437869 -1.03007054 3.87171217 xnew: 2.37644032 1.32837370 -1.03387071 3.87153902
## Iteration: 18 xold: 2.37644032 1.32837370 -1.03387071 3.87153902 xnew: 2.37632856 1.32398142 -1.03664119 3.87139375
## Iteration: 19 xold: 2.37632856 1.32398142 -1.03664119 3.87139375 xnew: 2.37624077 1.32076973 -1.03866215 3.87127761
## Iteration: 20 xold: 2.37624077 1.32076973 -1.03866215 3.87127761 xnew: 2.37617337 1.31842182 -1.04013699 3.87118742
## Iteration: 21 xold: 2.37617337 1.31842182 -1.04013699 3.87118742 xnew: 2.37612237 1.31670566 -1.04121363 3.87111869
## Iteration: 22 xold: 2.37612237 1.31670566 -1.04121363 3.87111869 xnew: 2.37608419 1.31545139 -1.04199976 3.87106696
## Iteration: 23 xold: 2.37608419 1.31545139 -1.04199976 3.87106696 xnew: 2.37605579 1.31453479 -1.04257386 3.87102835
## [1] 2.376 1.315 -1.043 3.871