i have 2 given matrices
a1 <- matrix(c(0.4092951, 0.1611806, 0.4283178, 0.001206529), nrow = 1) a2 <- matrix(c(0.394223557, 0.140443266, 0.463980790, 0.001352387), nrow = 1) i have initial matrix
b <- matrix(c(0.4095868, 0.1612955, 0.4286231, 0.0004946572, 0, 0.2732351, 0.7260891, 0.0006757670, 0, 0, 0.9909494, 0.0090505527, 0, 0, 0, 1), nrow = 4, byrow = t) i need update 'b' such that
a1 %*% b = a2 the above optimization problem objective function minimize
(a1 %*% b - a2) which drive value of sum(absolute value(a1 %*% b - a2)) zero, subject constraints:
lower triangle(b) = 0 ;
rowsum(b) = 1
## creating data vector a1 , a2 data = c(as.numeric(a1), as.numeric(a2)) ## objective function min_obj <- function(p){ ## creating matrix recreate 'b' p1 <- matrix(rep(0, 16), nrow = 4) k = 1 for(i in 1:nrow(p1)){ (j in 1:ncol(p1)){ if(j >= i){ p1[i,j] <- p[k] k = k+1 } } } actual <- matrix(data[1:(length(data)/2)], nrow = 1) pred <- matrix(data[(length(data)/ 2 + 1):length(data)], nrow = 1) s <- (actual %*% p1) - pred sum(abs(s)) } ## initializing initial values b taking non-zero values init <- b[b>0] opt <- optim(init, min_obj, control = list(trace = t), method = "l-bfgs-b", lower = rep(0, length(init)), upper = rep(1, length(init))) transformed_b <- matrix(rep(0, 16), nrow = 4) k = 1 for(i in 1:nrow(transformed_b)){ (j in 1:ncol(transformed_b)){ if(j >= i){ transformed_b[i,j] <- opt$par[k] k = k+1 } } } transformed_b the issue transformed_b rowsum of matrix not 1. highly appreciated.
"optim" right choice. since row sums have 1, there 6 parameters, not 10 in attempt. diagonal uniquely determined values strictly above diagonal.
a1 <- matrix(c(0.4092951, 0.1611806, 0.4283178, 0.001206529), nrow = 1) a2 <- matrix(c(0.394223557, 0.140443266, 0.463980790, 0.001352387), nrow = 1) b <- matrix(c(0.4095868, 0.1612955, 0.4286231, 0.0004946572, 0, 0.2732351, 0.7260891, 0.0006757670, 0, 0, 0.9909494, 0.0090505527, 0, 0, 0, 1), nrow = 4, byrow = t) #====================================================================== # build upper triangular matrix rowsums 1: b <- function(x) { x <- matrix(c(0,x[1:3],0,0,x[4:5],0,0,0,x[6],rep(0,4)),4,4,byrow=true) diag(x) <- 1-rowsums(x) return(x) } #---------------------------------------------------------------------- # function want minimize: f <- function(x) { return (sum((a1%*%b(x) - a2)^2)) } #---------------------------------------------------------------------- #optimization: opt <- optim( par = c(b[1,2:4],b[2,3:4],b[3,4]), fn = f, lower = rep(0,6), method = "l-bfgs-b" ) optb <- b(opt$par) result:
> optb [,1] [,2] [,3] [,4] [1,] 0.9631998 0.03680017 0.0000000 0.0000000000 [2,] 0.0000000 0.77820700 0.2217930 0.0000000000 [3,] 0.0000000 0.00000000 0.9998392 0.0001608464 [4,] 0.0000000 0.00000000 0.0000000 1.0000000000 > a1 %*% optb - a2 [,1] [,2] [,3] [,4] [1,] 9.411998e-06 5.07363e-05 1.684534e-05 -7.696464e-05 > rowsums(optb) [1] 1 1 1 1 i chose sum of squares instead of sum of absolute values, since differentiable. makes easier "optim" find minimum, guess.
Comments
Post a Comment