Main menu
You are here
Problem filling standardized & subdiagonal matrix
I am trying to fill either a standardized or a subdiagonal (m x m) matrix with values read from an R matrix. The matrix has 26 rows & 26 columns, so there should be (m*(m-1))/2 = 325 elements in the subdiagonal part. I note that the documentation says:
When ‘type’ is either ‘Lower’, ‘Sdiag’, ‘Symm’, or ‘Stand’, and the arguments to ‘free’, ‘values’, ‘labels’, ‘lbound’, or ‘ubound’ are vectors with enough elements to populate exactly one half of the matrix, then mxMatrix() populates the lower triangle of the matrix (and transposes the lower triangle if the matrix is symmetric).
But this is really the wrong number of matrix elements to fill a subdiagonal or standardized matrix (or a symmetric matrix for that matter). The condition should be when the number of elements is m(m-1)/2 (Sdiag or Stand) or m(m+1)/2 (Symmetric). Furthermore, exactly half the number of matrix elements is difficult to provide when m is odd :).
Accordingly, I seem to be getting an error both from OpenMx and R, and I think it is a bug.
> is.matrix(mplusvec)
[1] TRUE
> dim(mplusvec)
[1] 325 1
> mxMatrix("Sdiag",nrow=26,ncol=26,values=mplusvec)
Error in sys.call(-10) : not that many frames on the stack
In case there was some miscalculation about whether the diagonal should be included, I tried to read into a smaller matrix, but the fails somewhat more explicitly:
> mxMatrix("Sdiag",nrow=25,ncol=25,values=mplusvec)
Error: Upper triangle or diagonal of values matrix in subdiagonal matrix 'untitled8' is not all zeros!
> mxMatrix("Sdiag",nrow=26,ncol=26,values=mplusvec)
Error in sys.call(-10) : not that many frames on the stack
- Log in or register to post comments
Printer-friendly version
Comments
#1
This behavior is the result of two bugs and a disagreement.
First bug: Error in sys.call(-10): not enough frames on the stack. That was the result of a dumb mistake I made during refactoring. I've cleaned up the error handling code.
Second bug: The documentation was misleading. It now says,
When ‘type’ is ‘Lower’ or ‘Symm’, then the arguments to ‘free’, ‘values’, ‘labels’, ‘lbound’, or ‘ubound’ may be vectors of length N * (N + 1) / 2, where N is the number of rows and columns of the matrix. When ‘type’ is ‘Sdiag’ or ‘Stand’, then the arguments to ‘free’, ‘values’, ‘labels’, ‘lbound’, or ‘ubound’ may be vectors of length N * (N - 1) / 2.The disagreement: Currently, if the user specifies a vector of length N * (N - 1) / 2, then we convert that vector into the appropriate matrix. Or, if the user specifies an matrix of length N * (N - 1) / 2, then we try to shove those values into a matrix of size N x N. Subsequently, the error checking fails because the resulting matrix is not a sub-diagonal matrix. I see two options to resolve this matter:
(A) If the input is a matrix of the wrong dimensions as the MxMatrix object, then we throw an error. This is the behavior that I advocate. I remember asking for this behavior a long time ago, and I was talked out of it. Now I wish I had held my ground.
(B) If the input is a matrix with N * (N - 1) / 2 elements, then we treat this matrix as a vector and convert the vector into the appropriate subdiagonal or standardized matrix of dimensions N x N. I believe Mike Neale is advocating this behavior. Sounds very dangerous to me.
#2
I am coming around to your way of thinking, for a couple of reasons. One, you're always right :). Two, feeding a vector into a matrix is a bit dodgy because sometimes one would like to do so with byrow=TRUE but there isn't really any such mxMatrix option (and perhaps with good reason). For symmetric matrices, it is possible to use the xpnd() function from MCMCpack, which is effectively byrow=FALSE. The following is a byrow=TRUE version
xpndbyrow <- function (x, nrow = NULL)
{
dim(x) <- NULL
if (is.null(nrow))
nrow <- (-1 + sqrt(1 + 8 * length(x)))/2
output <- matrix(0, nrow, nrow)
k <- 0
for (i in 1:nrow) {
for (j in 1:i) {
k <- k+1
output[i,j] <- x[k]
output[j,i] <- x[k]
}
}
return(output)
}
This covers most eventualities of filling a symmetric matrix from a 'compact' vector such as might be generated by vech(), or might be output by a program such as classic Mx or Mplus.
Feeding Standardized matrices with compact, missing the diagonal, vectors is a pain. Better to insist on having the 1's in the vector, then using xpnd() or xpndbyrow() to fill out the matrix on the fly, methinks. Yes writing a more generic function which has an argument byrow=FALSE would be a lot better!
Filling lower triangular matrices from vectors might also be done on the way in rather than 'automagically'.
So, I would accept a modification if we have a suite of functions to expand matrices which are supplied in compact form, with suitable explicit arguments as to how we would like the vectorized input to be expanded. omxExpand() anyone?
#3
I didn't understand the beginning of your previous comment,
We do have an option for feeding a vector into a matrix using byrow=TRUE. What we don't have is an option for feeding a matrix of one size into a matrix of another size, using byrow=TRUE. I guess we could dictate that under these circumstances, the input matrix is transformed into a vector, and then byrow=TRUE can be applied to place the vector into the output matrix. But this seems convoluted in my opinion.
#4
I recall the matrix option being added: it makes sense when you have a matrix of the right size already, but not when you have a vector. So I think we should require matrices to be nrow x ncol in size.
Except: If they are have either 1 column or 1 row (i.e., is just a vector in disguise).
I can't see a problem applying as.vector to the values parameter after checking it has (nrow==1|ncol==1)
The chance of someone accidentally sending a 1D matrix of the right length is rather low.
#5
mxMatrix could be ultra-picky and only accept matrices of the right size, period. This would probably break some existing code. If the error message directed the user to a generic function for reformatting matrices, omxReshape() or something, then it would be a bit more bulletproof. The real danger is if mxMatrix() does something unexpected by the user. Still, omxReshape() would in a sense be shifting the problem away from mxMatrix() and into another function, creating a sort of function creep. Which I am in danger of becoming myself :).
#6
I'm closing this issue. It can be reopened later if necessary.