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.
The following MATLAB script constructs the Q matrix using a convex optimization framework:
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
Objective: Maximize the trace of the top-left nxn block of Q, i.e., maximize trace(Q_{1:n,1:n}).
Constraints:
ind_diag, ind_bl)xQ) = R_{ij}^{k}, where:ind_diag: A diagonal matrix defining the block structure for lag k.ind_bl: A block matrix with a single nonzero element at position
(i,j).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))
}
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)
)
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)
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?
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)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)
}
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 Fixupsfdiag <- 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