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
- Q must be positive semi-definite.
- The block diagonals of Q are constrained to match the elements of R. Specifically, for every lag k (from 0 to p), the constraints ensure that the trace of each block matches the corresponding block of R.
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:
- Q >= 0 (positive semidefinite constraint).
- trace(kron(
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).
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?
1