Home > Net >  How can I find number of iterations from this function in R
How can I find number of iterations from this function in R

Time:10-18

Maybe you can help me with this code. I got only 4 values, which are right and needed, but also I need to make a table with all iterations, till I get those right values, but I don't know how to get number of those iterations from this function I have: example provided here enter image description here

E <- matrix(c(1, 0, 0, 0,
              0, 1, 0, 0, 
              0, 0, 1, 0, 
              0, 0, 0, 1), byrow = TRUE, nrow = 4)

B <- matrix(c(1.2, 2, 3, 1.5), nrow = 4, byrow = TRUE)

D <- matrix(c(0.18, 0.57, 0.38, 0.42,
              0.57, 0.95, 0.70, 0.44,
              0.38, 0.70, 0.37, 0.18,
              0.42, 0.44, 0.18, 0.40), byrow = TRUE, nrow = 4)

#my matrix
A <- D   0.1*(8 3)*E
A

# Define a convenience function matching `numpy.inner`
inner <- function(a, b) as.numeric(as.numeric(a) %*% as.numeric(b))

conjugate_grad <- function(A, b) {    
  n <- dim(A)[1]
  x <- numeric(n)
  z <- b - A %*% x
  p <- z
  z_old <- inner(z, z)
  for (i in 1:60) {
    teta <- z_old / inner(p, (A %*% p))
    x <- x   teta * p
    z <- z - teta * A %*% p
    z_new <- inner(z, z)
    if (sqrt(z_new) < 0.001)
      break
    beta <- z_new / z_old
    p <- z   beta * p
    z_old <- z_new
  }
  return(x) 
}

conjugate_grad(A,B)


Thank you in advance!

CodePudding user response:

I think this modification of your code gives you what you want:

conjugate_grad <- function(A, b) {    
  n <- dim(A)[1]
  x <- numeric(n)
  history <- NULL   # Insert this line 
  z <- b - A %*% x
  p <- z
  z_old <- inner(z, z)
  for (i in 1:60) {
    teta <- z_old / inner(p, (A %*% p))
    x <- x   teta * p
    history <- cbind(history, x)   # Insert this line
    z <- z - teta * A %*% p
    z_new <- inner(z, z)
    if (sqrt(z_new) < 0.001)
      break
    beta <- z_new / z_old
    p <- z   beta * p
    z_old <- z_new
  }
  return(history)      # return history instead of x
}

conjugate_grad(A,B)
#           [,1]      [,2]       [,3]       [,4]
# [1,] 0.4326431 0.1288703 0.09609509 0.08061088
# [2,] 0.7210718 0.1695202 0.15988971 0.16900513
# [3,] 1.0816077 1.8689058 1.85211220 1.85311415
# [4,] 0.5408039 0.6284172 0.70703113 0.70548042

You have to accumulate the results as you compute them in an object, here called history.

  • Related