Search code examples
rmatrixrestriction

How to create a positive definite matrix with restrictions?


Consider that we want to create a covariance matrix of three random variables, v1, v2, v3. Then, the matrix has to satisfy the positive definiteness as well as other covariance restrictions such as positive variances and range of correlation.

Additionally, I want to impose a restriction: -var(v3) = cov(v2, v3) - cov(v1, v3).

For example, I need a matrix like below:

> tmp <- matrix(c(1, 0.6, 0.3,
                  0.6, 1, -0.2,
                  0.3, -0.2, 0.5), ncol = 3, byrow = T)
> is.positive.definite(tmp)
[1] TRUE

Note that var(v3) = 0.5, cov(v1, v3) = 0.3, cov(v2, v3) = -0.2 and thus -var(v3) = cov(v2, v3) - cov(v1, v3).

Note also that the variances (diagonal elements) are positive and the correlation coefficients such as 0.3/(1*0.5) are in [-1,1].

Of course, it is possible to make these matrices manually, but, it is so time-consuming because I have to make about a thousand of this kind of matrices.

Thus, how to make create a positive definite matrix with restrictions?


Solution

  • Here is a way with the CVXR package. It does not allow to impose on a variable matrix to be positive definite, only positive semidefinite. So I add a restriction on the determinant. One can only add a restriction on the log-determinant, it should be log_det(A) > -Inf but this is not possible so I took >= 0 (so det(A) >= 1). Also, one needs an objective. I choose to minimize the norm of A`. There is an infinity of solutions to your problem so one needs such additional constraints.

    library(CVXR)
    
    A <- Variable(3, 3, PSD = TRUE)
    constraints <- list(
      A[3,3] == A[1,3] - A[2,3], 
      log_det(A) >= 0
    )
    objective <- Minimize(norm(A))
    prob <- Problem(objective, constraints)
    sol <- psolve(prob, solver = "SCS")
    sol$getValue(A)