CODEWITHSUNDEEP
rmatlabmathematical-optimizationcvxr
Saïd Maanan
Saïd Maanan

Reputation: 697

R: Error Running Convex Optimization With CVXR

I am translating a MATLAB script that performs spectral factorization into R using the CVXR package for convex optimization. I have simplified the problem to focus only on constructing the Q matrix while maintaining the core constraints and objective.

MATLAB Script

The following MATLAB script constructs the Q matrix using a convex optimization framework:

The Constraints

function Q = compute_Q(n, p, Y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%input:
% n     = number of variables
% p     = order of the true AR model
% Y     = causal part of the "true" matrix polynomial $S^{-1}$ (IS)
%output:
% Q     = symmetric positive semidefinite matrix satisfying the constraints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

R = Y;
s = n * (p + 1);

cvx_begin quiet
    variable Q(s, s) symmetric;
    maximize trace(Q(1:n, 1:n));
    subject to
        Q == semidefinite(s);
        for k = 0:p
            ind_diag = diag(ones(p + 1 - k, 1), k);
            for i = 1:n
                for j = 1:n
                    ind_bl = zeros(n);
                    ind_bl(i, j) = 1;
                    trace(kron(ind_diag, ind_bl) * Q) == R{k + 1}(i, j);
                end
            end
        end
cvx_end

end

Problem

Objective: Maximize the trace of the top-left nxn block of Q, i.e., maximize trace(Q_{1:n,1:n}).

Constraints:

  1. Q >= 0 (positive semidefinite constraint).
  2. trace(kron(ind_diag, ind_bl)xQ) = R_{ij}^{k}, where:

R Translation

I translated the above MATLAB script into R using the CVXR package:

library(CVXR)

compute_Q <- function(n, p, Y) {
  R <- Y
  s <- n * (p + 1)
  Q <- Variable(s, s, symmetric = TRUE)
  
  # Objective: Maximize trace of the top-left block of Q
  objective <- Maximize(sum(diag(Q[1:n, 1:n])))
  
  # Constraints: Q is positive semidefinite
  constraints <- list(Q %>>% 0)
  
  # Add trace constraints for each block
  for (k in 0:p) {
    ind_diag <- diag(rep(1, p + 1 - k), k)
    for (i in 1:n) {
      for (j in 1:n) {
        ind_bl <- matrix(0, n, n)
        ind_bl[i, j] <- 1
        constraints <- c(
          constraints,
          trace(Q %*% kronecker(ind_diag, ind_bl)) == R[[k + 1]][i, j]
        )
      }
    }
  }
  
  # Solve the problem
  problem <- Problem(objective, constraints)
  result <- solve(problem, solver = "SCS")
  
  # Return the computed Q matrix
  return(result$getValue(Q))
}

Inputs

The inputs for the R function are:

k <- 9
n <- 5
p <- 1
Y <- list(
  matrix(c(2.1983699, 0.2066341, 0.2048484, 0.2066970, 0.2011900,
            0.2066341, 2.1954148, 0.2043585, 0.1901186, 0.1995567,
            0.2048484, 0.2043585, 2.2169782, 0.2075691, 0.2009464,
            0.2066970, 0.1901186, 0.2075691, 2.2021540, 0.2096051,
            0.2011900, 0.1995567, 0.2009464, 0.2096051, 2.1970326),
         nrow = 5, byrow = TRUE),
  matrix(c(0.1849041, 0.2060889, 0.1936582, 0.2057210, 0.2185530,
            0.1942851, 0.2074143, 0.1904800, 0.1798709, 0.2280653,
            0.1866000, 0.1944861, 0.1887652, 0.1946746, 0.1985805,
            0.2231732, 0.1974313, 0.2047517, 0.1884571, 0.1967431,
            0.2076642, 0.2123295, 0.1948372, 0.2092436, 0.1989898),
         nrow = 5, byrow = TRUE)
)

Error Encountered

When I run the R function:

Q_result <- compute_Q(n, p, Y)

I get the following error:

Error in .local(.Object, ...) : Invalid dimensions 00
17.
stop("Invalid dimensions ", dim)
16.
.local(.Object, ...)
15.
.nextMethod(.Object, ..., dim = intf_dim(.Object@value))
14.
eval(call, callEnv)
13.
eval(call, callEnv)
12.
callNextMethod(.Object, ..., dim = intf_dim(.Object@value))
11.
.local(.Object, ...)
10.
initialize(value, ...)
9.
initialize(value, ...)
8.
new(structure("Constant", package = "CVXR"), ...)
7.
.Constant(value = value)
6.
Constant(value = expr)
5.
as.Constant(y)
4.
Q %*% kronecker(ind_diag, ind_bl)
3.
Q %*% kronecker(ind_diag, ind_bl)
2.
trace(Q %*% kronecker(ind_diag, ind_bl))
1.
compute_Q(n, p, Y)

Question

Why does the CVXR implementation fail with this "Invalid dimensions" error? How should the constraints be set up in R to correctly mimic the MATLAB constraints?

Upvotes: 2

Views: 61

Answers (1)

ThomasIsCoding
ThomasIsCoding

Reputation: 102529

Reasons

There are two places you need to change in your R code, since their definitions are absolutely different from the Matlab functions

  • diag: In R, the argument k in diag(v,k) refers to the number of rows of the diagonal matrix, instead of the expected "the position of diagonal" (as done by Matlab)
  • trace: In R, the function trace is a call to interactively keep track of debugging, rather than computing the trace of a matrix (as done by Matlab)

Fixups

  • regarding diag, you should define a customized diag function to mimic the version in Matlab, e.g.,
fdiag <- function(v, k) {
  n <- length(v) + abs(k)
  mat <- matrix(0, n, n)
  replace(mat, col(mat) - row(mat) == k, v)
}
  • Regarding trace, you should replace it by matrix_trace, which is from CVXR package itself and is dedicated to compute the trace of a function, as desired.

R Code with Above Fixups

fdiag <- function(v, k) {
  n <- length(v) + abs(k)
  mat <- matrix(0, n, n)
  replace(mat, col(mat) - row(mat) == k, v)
}


library(CVXR)
compute_Q <- function(n, p, Y) {
  R <- Y
  s <- n * (p + 1)
  Q <- Variable(s, s, symmetric = TRUE)

  # Objective: Maximize trace of the top-left block of Q
  objective <- Maximize(sum(diag(Q[1:n, 1:n])))

  # Constraints: Q is positive semidefinite
  constraints <- list(Q %>>% 0)

  # Add trace constraints for each block
  for (k in 0:p) {
    ind_diag <- fdiag(rep(1, p + 1 - k), k)
    for (i in 1:n) {
      for (j in 1:n) {
        ind_bl <- matrix(0, n, n)
        ind_bl[i, j] <- 1
        constraints <- c(
          constraints,
          matrix_trace(Q %*% kronecker(ind_diag, ind_bl)) == R[[k + 1]][i, j]
        )
      }
    }
  }

  # Solve the problem
  problem <- Problem(objective, constraints)
  result <- solve(problem, solver = "SCS")

  # Return the computed Q matrix
  return(result$getValue(Q))
}

gives

> (Q <- compute_Q(n, p, Y))
           [,1]      [,2]      [,3]      [,4]      [,5]       [,6]       [,7]
 [1,] 2.1225066 0.1298917 0.1313423 0.1330286 0.1229291 0.18490392 0.20608872
 [2,] 0.1298917 2.1167175 0.1294767 0.1147492 0.1192947 0.19428491 0.20741411
 [3,] 0.1313423 0.1294767 2.1451438 0.1356292 0.1244219 0.18659984 0.19448594
 [4,] 0.1330286 0.1147492 0.1356292 2.1293226 0.1327116 0.22317303 0.19743112
 [5,] 0.1229291 0.1192947 0.1244219 0.1327116 2.1144291 0.20766400 0.21232930
 [6,] 0.1849039 0.1942849 0.1865998 0.2231730 0.2076640 0.07586326 0.07674252
 [7,] 0.2060887 0.2074141 0.1944859 0.1974311 0.2123293 0.07674252 0.07869718
 [8,] 0.1936580 0.1904798 0.1887651 0.2047515 0.1948370 0.07350616 0.07488188
 [9,] 0.2057208 0.1798707 0.1946745 0.1884569 0.2092434 0.07366851 0.07536945
[10,] 0.2185528 0.2280651 0.1985803 0.1967429 0.1989896 0.07826099 0.08026212
            [,8]       [,9]      [,10]
 [1,] 0.19365803 0.20572083 0.21855281
 [2,] 0.19047982 0.17987071 0.22806510
 [3,] 0.18876506 0.19467445 0.19858034
 [4,] 0.20475154 0.18845692 0.19674291
 [5,] 0.19483701 0.20924341 0.19898959
 [6,] 0.07350616 0.07366851 0.07826099
 [7,] 0.07488188 0.07536945 0.08026212
 [8,] 0.07183425 0.07193995 0.07652456
 [9,] 0.07193995 0.07283130 0.07689358
[10,] 0.07652456 0.07689358 0.08260348

Upvotes: 2

Related Questions