implement [] in mxAlgebra/mxConstraint

tbates's picture
Project:OpenMx
Component:Code
Category:feature request
Priority:normal
Assigned:Unassigned
Status:closed
Description

[ (the bracket accessor function) is not listed as an implemented function in mxAlgebra.

It's handy, but not available in algebra or the algebra of a constraint:

# ======================================
# = 1. Try bracket access in mxAlgebra =
# ======================================
nOrd = 3
nVar = 5
fit1 = mxModel("fit1",
	mxMatrix(name="a"   , "Full", nrow=nVar, ncol=nVar, free=T, values=1:25),
	mxMatrix(name="unit", "Unit", nrow=nOrd, ncol=1),
	mxAlgebra(diag2vec(a)[1:3], name = "ordElements"),
	mxConstraint(name = "con1", ordElements == unit),
	mxFitFunctionAlgebra("ordElements")
)
fit1 = mxRun(fit1)
# Error: In algebra 'fit1.ordElements' could not find function with name [ and 2 arguments.
 
# =========================================
# = 1. Try bracket access in mxConstraint =
# =========================================
fit1 = mxModel("fit1",
	mxMatrix(name="a"   , "Full", nrow=nVar, ncol=nVar, free=T, values=1:25),
	mxMatrix(name="unit", "Unit", nrow=nOrd, ncol=1),
	mxAlgebra(diag2vec(a), name = "diagElements"),
	mxConstraint(name = "con1", diagElements[1:3] == unit),
	mxFitFunctionAlgebra("ordElements")
)
fit1 = mxRun(fit1)
# Error: In algebra 'fit1.untitled52' could not find function with name [ and 2 arguments.

One kind of work-around is that Bracket access of a single cell "works" because it is interpreted as a bracket-style matrix address, so this is OK:

# Goal = constrain third cell on diagonal of "a" == 1
fit1 = mxModel("fit1",
	mxMatrix(name="a", "Full", nrow=5, ncol=5, free=T, values=1:25),
 	# next line fails with error
	mxConstraint(a[1,1] == 1, name = "out1"),
	mxConstraint(a[2,2] == 1, name = "out2"),
	mxConstraint(a[3,3] == 1, name = "out3"),
 	mxFitFunctionAlgebra("a")
)
fit1 = mxRun(fit1); mxEval(a, fit1)
 
# Running fit1 
# Warning message:
# In model 'fit1' Optimizer returned a non-zero status code 1. The final iterate satisfies the optimality 
# conditions to the accuracy requested, but the sequence of iterates has not yet converged. Optimizer was
# terminated because no further improvement could be made in the merit function (Mx status GREEN). 
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    1   12   17   22
[3,]    3    8    1   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25

There's a work-around of using matrix transforms to do the accessing, but easier to avoid if possible? ...

# ====================================
# = Use matrix transform work around =
# ====================================
 
nOrd = 3
nVar = 5
fit1 = mxModel("fit1",
	mxMatrix(name="a", "Full", nrow=nVar, ncol=nVar, free=T, values=1:25),
	mxMatrix(name="b", "Full", nrow=nVar, ncol=nOrd, free=T, values=
		c(1, 0, 0,
		  0, 1, 0,
		  0, 0, 1,
		  0, 0, 0,
		  0, 0, 0
		)
	),
	mxMatrix(name="unit", "Full", nrow=1, ncol=nOrd, free=F, values=1),
	mxAlgebra(t(diag2vec(a)) %*% b, name="z"),
	mxConstraint(name = "f", z == unit),
	mxFitFunctionAlgebra("z", numObs= NA, numStats=NA)
)
fit1 = mxRun(fit1)
# Running fit1
# Warning message:
# In model 'fit1' Optimizer returned a non-zero status code 1. The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. Optimizer was terminated because no further improvement could be made in the merit function (Mx status GREEN).
mxEval(list(z), fit1)
[[1]]
          [,1]      [,2] [,3]
[1,] 0.9999998 0.9999998    1

Comments

mhunter's picture

#1

This seems very difficult to implement, and there are several ways around it. My favored method is your "matrix transform". I just call it filtering, and I use the handy diag() function in R to get the nice values. Here's my version.

nOrd = 3
nVar = 5
fit1 = mxModel("fit1",
	mxMatrix(name="a"   , "Full", nrow=nVar, ncol=nVar, free=T, values=1:25),
	mxMatrix(name="unit", "Unit", nrow=nOrd, ncol=1),
	mxMatrix(name="Filter", "Full", nrow=nOrd, ncol=nVar, free=FALSE, values=diag(1, nrow=nOrd, ncol=nVar)),
	mxAlgebra(Filter %*% diag2vec(a), name="ordElements"),
	mxConstraint(name = "con1", ordElements == unit),
	mxFitFunctionAlgebra("ordElements")
)
fit1 = mxRun(fit1)
mxEval(list(ordElements), fit1)
#> fit1 = mxRun(fit1)
#Running fit1 
#> mxEval(list(ordElements), fit1)
#[[1]]
#     [,1]
#[1,]    1
#[2,]    1
#[3,]    1

tbates's picture

#2

Yip, can use filtering as a workaround (diag shortcut won't work for the more arbitrary arrangements, but will write a helper to do this work in one line.

Just imagining this in R though: "Oh you want to get and set object elements? yeah, you have to create custom matrix algebra for that"... :-) :-(

mhunter's picture

#3

Yeah, I can see your point that this seems elementary. But algebras are a relatively advanced feature. Arbitrary (i.e. fluent replication of R behavior) square-bracket usage within algebras is even more advanced. And it's more "Oh you want to get and set object elements? Here are six ways you can do it, and one way you can't.".

tbates's picture

#4

It's interesting how very quickly one requires advanced features in OpenMx. Perhaps that's part of the issue.

For instance it's hard to make a research-ready twin script that doesn't require mxAlgebras, mxConstraints, definition variables, and, in most "real data" threshold matrix construction, constraining of algebras stepping around a mixture of binary, ordinal and continuous variables, plus fairly sophisticated start-value regimes to get the optimizer off the ground...

sic itur ad astra :-)

mhunter's picture

#5

I don't view twin modeling as beginner-level. To me this is a reason why BGA is filled with so many method-savvy researchers: they have to be methodologically sophisticated to be doing BG. I don't think these models are easier to specify in another language. This reinforces my thinking that it's not that OpenMx makes these simple models hard, rather it's that these are hard models. I grant you that we should endeavor to make model specification easier. In any case, it sounds like a nice opportunity for helper functions.

One of the things the Dev Team has discussed is a family of model generation helper functions for various often-used models: factor model, latent growth curve, ACE, etc. This could be a feature request managed among the Issues here. You could add a request for a helper function for the kind of model you're discussing.

mhunter's picture

#6

Status:active» closed

These things don't work By Design.