another mxConstraint question

8 replies [Last post]
manu's picture
Offline
Joined: 10/08/2009

Dear all,

I have a large matrix (A), and I want to constrain certain submatrices within A to (let's say the product of) some smaller matrices. How can this be done in OpenMx?

I am looking for something like:

A <- mxMatrix("Full", 4, 4, labels = c("a","b",NA,NA,"c","d",NA,NA,NA,NA,"a","b",NA,NA,"c","d"), byrow = TRUE, free = TRUE, name = "AA")

B <- mxMatrix("Full", 2, 2, labels = c("w","x","y","z"), byrow = TRUE, free = TRUE, name = "BB")

C <- mxMatrix("Full", 2, 2, labels = c("w1","x1","y1","z1"), byrow = TRUE, free = TRUE, name = "CC")

mxConstraint(A[1:2,1:2] == B%*%C)
mxConstraint(A[3:4,3:4] == B%*%C)

OR

mxConstraint(AA[1:2,1:2] == BB%*%CC)
mxConstraint(AA[3:4,3:4] == BB%*%CC)

Part of the problem appears to be that I am mixing up R and openMx Code.

Thank you very much for your help!

mspiegel's picture
Offline
Joined: 07/31/2009
The second form should work

The second form should work in your model:

mxConstraint(AA[1:2,1:2] == BB %*% CC, "constraint1")
mxConstraint(AA[3:4,3:4] == BB %*% CC, "constraint2")

If the names assigned to the matrices were "A", "B", and "C", then you could use the form:

mxConstraint(A[1:2,1:2] == B %*% C, "constraint1")
mxConstraint(A[3:4,3:4] == B %*% C, "constraint2")

One small optimization would be to declare an algebra that yielded B %*% C in your model, and then use that algebra on the right-hand side of both constraints. OpenMx does not yet perform common subexpression elimination.

manu's picture
Offline
Joined: 10/08/2009
thank you!!! I still do not

thank you!!!
I still do not get it to run but I keep on eliminating errors;-)
One of the "errors" I eliminated in my skript was the following: In R scalar multiplication is simple. E.g.:
a <- matrix(c(1,2,3,4), nrow=2, ncol=2)
2*a
[,1] [,2]
[1,] 2 6
[2,] 4 8

This is different in OpenMx:

test <- mxModel("test",
mxMatrix(values=a, name = "a"),
mxAlgebra(2*a, "result"))
out <- mxRun(test)
out@algebras$result

does not do the scalar multiplication (but also does not return an error message), so one has to define a matrix first and do element-multiplication:

test <- mxModel("test",
mxMatrix(values=a, name = "a"),
mxMatrix(type="Full", nrow=2, ncol=2, values =c(2,2,2,2), name = "b"),
mxAlgebra(b*a, "result"))
out <- mxRun(test)
out@algebras$result

It took me a while to figure that out, so I thought I share it with other users, who may stumble over the same problem.

mspiegel's picture
Offline
Joined: 07/31/2009
There is a confluence of

There is a confluence of events in your example. The '*' operator executes elementwise multiplication, as you have discovered. It should be that the first model throws an error because the two matrices are not conformable. We changed the mechanism for issuing errors, and apparently we introduced a bug in our new mechanism. This will be fixed in the next binary release. You can see the error message by looking at out@output$status[[3]]. But it should have thrown an error, sorry about that.

Also, note there are three types of multiplication in OpenMx. They are dot product (element-wise multiply), matrix multiplication, and kronecker product. See http://openmx.psyc.virginia.edu/wiki/matrix-operators-and-functions for a table. The following script will do the right thing:

test <- mxModel("test",
mxMatrix(values=a, name = "a"),
mxAlgebra(2 %x% a, "result"))
out <- mxRun(test)

UPDATE: breath easy, almost all errors are still detected in the current OpenMx binary release. We execute all the algebras in R, and catch 90% of the errors before we start executing the optimizer. The reason this error was not detected is that the '*' operator in R is used for both element-wise matrix-matrix multiplication and for scalar-matrix multiplication. In the back-end, all operations are matrix operations. There do not exist scalar or vector types in the OpenMx backend. Hence this specific error must be detected in the backend.

manu's picture
Offline
Joined: 10/08/2009
Thank you - after some forth

Thank you - after some forth and back I figured that would be the case and is easy to resolve.
Unfortunately my model still does not run and I think the problem is related to Jan's post a few days ago (http://openmx.psyc.virginia.edu/thread/519).
Assume I formulate a model where (I-A)^-1 is singular, thus the model cannot be estimated. Only after introducing a few very simple equality constraints of the type
a == 0.5*b
c ==1 + 0.5*b
the model becomes nonsingular and runs smoothly in Mplus. When imposing the constraints via mxConstraint in OpenMx, however, I continue to get the error message: " The algebra 'Z' in model 'xyz' generated the error message: Lapackroutine dgesv: system is singular". So it appears that OpenMx throws out the error message before taking my equality constraints into account. Is that true? In the manual it is recommended that the mxConstraint function should not be used to constrain free parameters, but that one should assign the same labels. This, however, appears to be no solution in my case where a and b cannot be assigned the same label since a == 0.5*b. Am I missing something? Do you have any recommendations? I know, this is not a very clear-cut question, but I am trying to understand how mxConstraint() works and what goes wrong in my case. Once again, thank you very much for your help!

PS. Ideally I would perfer to "equate" submatrices instead of single parameters as noted in my first post, but this appears to be unrelated to the actual problem.

tbrick's picture
Offline
Joined: 07/31/2009
If I were to guess why you

If I were to guess why you were getting this error, I would guess that your starting values for this model don't meet your constraints or that they cause your matrix to be singular for other reasons.

This often happens if you're inverting a matrix of free parameters, and you set them all to identical starting values. A statement like:
 mxMatrix("Full", 2, 2, free=T, values = rep(.8, 4) + rnorm(4, sd=.1))
will give you starting values with random offsets, and make your matrix invertible.

There is also a way to use labels to enforce these constraints. Consider:

 mxAlgebra(.5 * b, name="a")
mxAlgebra(1 + .5 * b, name="c")

Now you can reference "a[1,1]" as a label the same way you would reference the free parameter a in your original model.

This approach will also extend to matrices or submatrices (you can access pieces using the : and [] operators).

manu's picture
Offline
Joined: 10/08/2009
Thank you - it finally

Thank you - it finally works:-)
As you said, the starting values caused the model (I-A) to be singular. I knew that, but did not care because I constrained the critical parts of the A matrix to another parameterized matrix with correct starting values and assumed the initial values would simply be overwritten by mxConstraint. This, however, is not the case, so I kept getting the nonsingular error message. The solution was to set the relevant starting values to NA and free=FALSE, before constraining them to another matrix with correct starting values and free=TRUE. In hindsight this makes perfect sense...
Using mxAlgebra and referencing the labels as you described above is perfect for that purpose - this is a wonderful option!

mspiegel's picture
Offline
Joined: 07/31/2009
"The solution was to set the

"The solution was to set the relevant starting values to NA and free=FALSE, before constraining them to another matrix with correct starting values and free=TRUE."

I just want to clarify whether you're using a mxConstraint() to create an equality constraint, or using a mxAlgebra() and referencing the labels? Because I don't think using an mxConstraint() with starting values of NA will work. Check the final values of those cells after optimization, and confirm they are not NA values.

A little bit of background: a numeric NA value in R is translated into a IEEE floating point value of NaN when passed to a function in C. And the NaN value has the property that (NaN == x) returns FALSE for all values of x, including the value of NaN. Which is why an equality constraint involving a NA is suspicious.

manu's picture
Offline
Joined: 10/08/2009
yes, I use mxAlgebra() and

yes, I use mxAlgebra() and reference the labels and that gives me correct results. After setting the starting values to NA I did not try it with mxConstraint() -- and I am not going to try it now, since your explanation makes perfect sense. Thank you!