CVXR Elementwise multiplication of matrix with vector

336 views Asked by At

I want to use CVXR to find the optimal values of a vector. In the objective function i need to multiply a matrix with a vector in an elementwise way:

b: Nx1 vector X: Nxp matrix result: Nxp matrix

Example:

# Set the dims
N <- 50
p <- 5
X <- matrix(rnorm(N*p), N, p)

# to find the optimal values using optim() one could simply have a numeric object
# say the optimal values are 0.1, -0.2, 0.3, -0.5, 0.6
b <- c(0.1, -0.2, 0.3, -0.5, 0.6)

# Then we can have the Nxp matrix with the product
# (where column i of X is multiplied by element i of b) is given by
X*b

b is the vector of coefficient to be optimised.

Using CVXR one must declare

b <- Variable(p)

as Variable object with uses the matrix form, therefore later we can't really multiply as in the previous case.

Also, we don't want to create a matrix of b: Nxp because we need to have one optimal value for all N observations of the i-th column (therefore the mul_elemwise(X, X*b) option wouldn't work as it would give different optimal values for different observations of N - if I am not mistaken).

thanks,

1

There are 1 answers

0
Erwin Kalvelagen On

To recap: this is the R behavior:

> m <- 2
> n <- 3
> A <- matrix(c(1,2,3,4,5,6),m,n)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> x <- c(1,2)
> A*x
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    4    8   12
> 

This is essentially

A[i,j]*x[i]

R behind the scenes extends (recycles) x to have it as many elements as A and then does elementwise multiplication in a columnwise fashion.

In CVXR things are a bit different. %*% is for matrix multiplication and * is for elementwise multiplication. But CVXR does not do this recycling. So for A*x it requires A and x to have the same shape (i.e. an (mxn) matrix).

This means we need to do this extending (recycling) ourselves. This can be written as follows:

> x %*% t(rep(1,n))
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2

So we can write:

> A * (x %*% t(rep(1,n)))
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    4    8   12

This is what we can use in the CVXR model:

> library(CVXR)
> x <- Variable(m)
> Y <- Variable(m,n)
> e <- t(rep(1,n))
> e
     [,1] [,2] [,3]
[1,]    1    1    1
> problem <- Problem(Minimize(0),list(x == c(1,2), Y == A * (x %*% e)) )
> sol <- solve(problem)
> sol$status
[1] "optimal"
> sol$getValue(x)
     [,1]
[1,]    1
[2,]    2
> sol$getValue(Y)
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    4    8   12
>