Using Conjugate Gradient to Solve Matrix Equations with R-CodePudding

I would really appreciate your help with Conjugate Gradient method to Solve Matrix Equations with R. I found similar solution with python, but I don't know how to write exact code with R. enter image description here

import numpy as np

def conjugate_grad(A, b, maxiter = 5):
    n = A.shape[0]
    x = np.zeros(n)
    r = b - A @ x
    p = r.copy()
    r_old = np.inner(r, r)
    for it in range(maxiter):
        alpha = r_old / np.inner(p, A @ p)
        x  = alpha * p
        r -= alpha * A @ p
        r_new = np.inner(r, r)
        if np.sqrt(r_new) < 1e-10:
            break
        beta = r_new / r_old
        p = r   beta * p
        r_old = r_new.copy()
    return x

import numpy as np

A = np.array([[4, 1], [1, 3]])
b = np.array([1, 2])
x = conjugate_grad(A, b)

The algorithm of Conjugate Gradient method is like thisenter image description here

I have started writing code in R, but I know that something is wrong with it,I don't know how to continue after break.

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



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

exactness <- NULL

epsilon <- 0.0001



if (isSymmetric(A) == TRUE) {
  print(paste("Matrix is symmetric"))
} else {print ("matrix is not symmetric")}


library(matrixcalc)
if (is.positive.definite(A) == TRUE) {
  print(paste("matrix is  positive-definite "))
} else {print ("Matrix is not  positive-definite")}

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

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


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


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

scalarmult <- NULL
scalarmult2 <- NULL
scalarmult3 <- NULL
scalarmult4 <- NULL
teta <- NULL
beta<-NULL


for (i in 1:60) {
  p[, i 1] <- A %*%x[,i] - B
  z[, i 1] <- A %*%x[,i] - B
  scalarmult[i 1] <- t(z[, i 1])%*%p[,i 1]
  r[, i 1] <- A%*%p[,i 1]
  scalarmult2[i 1] <- t(r[, i 1])%*%p[,i 1]
  teta[i 1] <- scalarmult[i 1]/scalarmult2[i 1]
  x[,i 1] <- x[,i]-teta[i 1]%*%t(p[,i 1])
  
  if (scalarmult[i 1] < epsilon^2) {
    print(paste("iteration number", i, ", exactness: ", scalarmult[i 1]))
    print(paste(x[,i 1]))
    break
    beta[i 1]<-
    
    
  }
}

Thank You in advance!

CodePudding user response：

I struggle to understand your code attempt. A one-to-one, literal translation of the Python code into R would be as follows:

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

conjugate_grad <- function(A, b, maxiter = 5) {    
    n <- dim(A)[1]
    x <- numeric(n)
    r <- b - A %*% x
    p <- r
    r_old <- inner(r, r)
    for (it in seq_len(maxiter)) {
        alpha <- r_old / inner(p, (A %*% p))
        x <- x   alpha * p
        r <- r - alpha * A %*% p
        r_new <- inner(r, r)
        if (sqrt(r_new) < 1e-10)
            break
        beta <- r_new / r_old
        p <- r   beta * p
        r_old <- r_new
    }
    return(x) 
}

Using the same sample data as in the Python example gives

A <- matrix(c(4, 1, 1, 3), ncol = 2)
b <- c(1, 2)
conjugate_grad(A, b)
#           [,1]
#[1,] 0.09090909
#[2,] 0.63636364

which matches the output of the Python routine and R's solve

solve(A, b)
#[1] 0.09090909 0.63636364

Please note that this a literal translation of the Python routine, which is rarely a good idea (as it's most likely not the most efficient and R-optimal solution). I leave it up to you to optimise this.