Error of using vech() in mxAlgebra()

4 replies [Last post]
Mike Cheung's picture
Offline
Joined: 10/08/2009

Hi, all.

Thanks for creating OpenMx, a wonderful R package!

According to the manual of version 0.2.3-1006, vech() is supported in mxAlgebra() via MCMCpack while vech() is not listed in http://openmx.psyc.virginia.edu/wiki/matrix-operators-and-functions . I tried it in version 0.2.3-1006 but it failed. This also applied to other functions, such as diag(), c() and c(t()) (a ")" was missing in the manual).

Here are the example.
> library(OpenMx)
> library(MCMCpack)
> mxVersion()
[1] "0.2.3-1006"
> A <- mxMatrix("Symm", nrow = 3, ncol = 3, values=1:6, name = "A")
> B <- mxAlgebra(vech(A), name = "B")
> model <- mxModel("test", A, B)
> fit <- mxRun(model)
Running test
Error in substituteOperators(as.list(retval)) :
Could not find function with name vech and 1 arguments

By the way, would it also be possible to implement a similar function of vech() that does not include the diagonals? This function would be useful for analyzing correlation structure. Thanks.

Regards,
Mike

mspiegel's picture
Offline
Joined: 07/31/2009
Oops. You've found a mistake

Oops. You've found a mistake in the mxAlgebra() documentation. The wiki list of operators and functions is up to date. I've checked in a patch to the source code repository; the correction will appear in the next binary release. I could imagine vech() with and without diagonals could be useful. Is it possible to generate the same functionality using the MxMatrix types 'Lower' and 'Sdiag'? Maybe perform a matrix multiplication on the symmetric matrix to extract the lower triangle (with or without diagonal)? In any case, thanks for alerting us to the mistake.

Mike Cheung's picture
Offline
Joined: 10/08/2009
Thanks for the suggestions.

Thanks for the suggestions. Since the vectorized matrix can be a computed matrix (e.g., an implied correlation/covariance matrix), I am not sure if the MxMatrix types 'Lower' and 'Sdiag' can be used to vectorize the matrix here.

Anyway, I have implemented an ugly version of it.
# Half-vectorization (vech) by column of a MxMatrix object
# X: MxMatrix object
# name: name of the half-vectorized matrix
vech_mx <- function(X, name) {
expr <- NULL
for (j in 1:ncol(X)) {
for (i in j:nrow(X)) {
expr <- paste(expr, X@name, "[", i, ",", j, "],", sep="")
}
}
expr <- sub(",$", ")", expr) # replace the last "," with ")"
expr <- paste(name, " <- mxAlgebra(rbind(", expr, ", name=\"", name, "\")", sep="")
eval(parse(text=expr))
}

# Example
A <- mxMatrix("Full", ncol=3, nrow=3, values=1:9, name="A")
vech_A <- vech_mx(A, "vech_A")
model <- mxModel("Test", A, vech_A)
fit <- mxRun(model)
mxEval(vech_A, fit)

mspiegel's picture
Offline
Joined: 07/31/2009
The functions vech() and

The functions vech() and vechs() have been implemented in the OpenMx library. They support half-vectorization and strict half-vectorization, respectively. I'll add documentation for both of these functions to the source code repository tomorrow. The functions will be available in the next binary release. To access them before the next release, use the HOWTO build OpenMx from source.

Mike Cheung's picture
Offline
Joined: 10/08/2009
It works. Many thanks!

It works. Many thanks!