Constrained Optimization in R not giving expected results

Question

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 an 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 to update 'b' such that

    a1 %*% b = a2

The above is an optimization problem where the objective function is to minimize

   (a1 %*% b - a2)

which would drive the value of the sum(absolute value(a1 %*% b - a2)) to zero, subject to the constraints:
Lower triangle(b) = 0 ;
RowSum(b) = 1

   ## creating a data vector with a1 and a2 
   data = c(as.numeric(a1), as.numeric(a2))

   ## objective function
   min_obj <- function(p){
   ## Creating a matrix to recreate 'b'
      p1 <- matrix(rep(0, 16), nrow = 4)
      k = 1
      for(i in 1:nrow(p1)){
          for (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 the initial values for b taking only 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)){
       for (j in 1:ncol(transformed_b)){
           if(j >= i){
              transformed_b[i,j] <- opt$par[k]
             k = k+1
           }
       }
    }

    transformed_b

The issue with transformed_b is that rowSum of the matrix is not 1. Any help is highly appreciated.

It would be helpful to also state what value you're receiving for the rowSum instead. — cryptic_star

mra68 mra68 · Accepted Answer · 2015-09-25T22:23:36

"optim" is the right choice. Since the row sums have to be 1, there are only 6 parameters, not 10 as in your attempt. The diagonal is uniquely determined by the values strictly above the 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 an upper triangular matrix with 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)
}

#----------------------------------------------------------------------
# The function we want to 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 the sum of squares instead of sum of absolute values, since it is differentiable. This makes it easier for "optim" to find the minimum, I guess.

Constrained Optimization in R not giving expected results

1 Answers