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
.