Classic Mx had a \chol function to obtain L for symmetric positive (semi-)definite A in A = L * L'. I figure that the BLAS function
dpbtrf could be used: SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
and that it wouldn't be too hard to add this to the list of omxAlgebraFunctions.c - but that it would be harder for me than for say Tim Brick to do this. I made a start but have been hauled off to do other frantically urgent things. I suspect we need to tell the algebra cruncher routine the priority of this operation also.
void omxMatrixCholesky(omxMatrix** matList, int numArgs, omxMatrix* result)
{
if(OMX_DEBUG_ALGEBRA) { Rprintf("ALGEBRA: Matrix Cholesky Decomposition.\n");}
omxMatrix* inMat = matList[0];
int ipiv[inMat->rows], lwork;
lwork = 4 * inMat->rows * inMat->cols;
double work[lwork];
int l = 0;
if(rows != cols) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Cholesky decomposition of non-square matrix cannot even be attempted\n");
omxRaiseError(result->currentState, -1, errstr);
Free(errstr);
}
omxCopyMatrix(result, inMat);
F77_CALL(dpbtrf)('L', &(result->cols), result->data, &(result->leading), ipiv, &l);
if(l != 0) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Attempted to Choleskify non-invertable matrix.");
omxRaiseError(result->currentState, -1, errstr);
Free(errstr);
}
}
I did the following:
cholesky <- function(x) {
x <- as.matrix(x)
if(nrow(x) != ncol(x)) {
stop("matrix must be square")
}
}
57 "Cholesky Decomposition" "cholesky" 1 "omxCholesky"
void omxCholesky(omxMatrix** matList, int numArgs, omxMatrix* result)
{
}
(which may or may not work). Everything compiles ok, but when I run the following code:
library(OpenMx)
C<-mxModel("C",A,B)
A<-mxMatrix('Symm',nrow=2,ncol=2,values=c(1,.5,.5,1),name="A")
B<-mxAlgebra(cholesky(A),name="algB")
C<-mxModel("C",A,B)
mxRun(C)
It complains with:
Error: The following error occurred while evaluating the subexpression 'cholesky(C.A)' during the evaluation of 'algB' in model 'C' : could not find function "cholesky"
So I feel like I must have missed something in connecting up the dots...
On second thought, is it true that you want the behavior of
cholesky(A)
to be identical to calling the R functionchol(A)
?If that's the case, then don't define the function
cholesky(A)
in R. Just use the R functionchol(A)
, and then define the function omxCholesky in C so that the back-end knows what to do. I would change your directions to:57 "Cholesky Decomposition" "chol" 1 "omxCholesky"
And then you don't need to change the NAMESPACE file. Besides using an existing R function to provide OpenMx functionality, this also helps in writing test cases, as the following should be identical:
> library(OpenMx)
>
> A<-mxMatrix('Symm',nrow=2,ncol=2,values=c(1,.5,.5,1),name="A")
> B<-mxAlgebra(chol(A),name="algB")
> C<-mxModel("C",A,B)
> runC <- mxRun(C)
Running C
> runC$algB
mxAlgebra 'algB'
@formula: chol(A)
@result:
[,1] [,2]
[1,] 1 0.5000000
[2,] 0 0.8660254
dimnames: NULL
>
> # test 2, 3x3
> A<-mxMatrix('Symm',nrow=3,ncol=3,values=c(1,.5,.5,.5,1,.5,.5,.5,1),name="A")
> B<-mxAlgebra(chol(A),name="algB")
> C<-mxModel("C",A,B)
> runC <- mxRun(C)
Running C
> runC$algB
mxAlgebra 'algB'
@formula: chol(A)
@result:
[,1] [,2] [,3]
[1,] 1 0.5000000 0.5000000
[2,] 0 0.8660254 0.2886751
[3,] 0 0.0000000 0.8164966
dimnames: NULL
>
> # Alternatively
> A<-mxMatrix('Symm',nrow=3,ncol=3,values=c(1,.5,.5,.5,1,.5,.5,.5,1),name="A")
> B <- mxAlgebra(chol(A), name = 'B')
> model <- mxModel('model', A, B)
> frontendCompute <- mxEval(B, model, compute = TRUE) # calculation in R
> backendCompute <- mxEval(B, mxRun(model), compute = FALSE) # calculation in backend
Running model
> frontendCompute
[,1] [,2] [,3]
[1,] 1 0.5000000 0.5000000
[2,] 0 0.8660254 0.2886751
[3,] 0 0.0000000 0.8164966
> backendCompute
[,1] [,2] [,3]
[1,] 1 0.5000000 0.5000000
[2,] 0 0.8660254 0.2886751
[3,] 0 0.0000000 0.8164966
I would prefer that chol(A) returns the lower triangle (consistent with classic Mx) but can deal with t(chol(A)) to get the desired result.
Should I upload these changes? I get a make test pass everything except for:
Number of errors: 1
From model models/passing/IntroSEM-ThreeLatentMultipleRegTest2.R :
[1] "In omxCheckCloseEnough(expectSE, as.vector(threeLatentMultipleReg1Out@output$standardError), 0.001) : '0.0769 0.087865 0.102543 0.127828 0.088084 0.119791 0.111236 0.152625 0.115206 0.113084 0.110164 0.122883 0.118535 0.114589 0.145932 0.10598 0.106799 0.141507 0.133841 0.171978 0.13313 0.109413 0.121718 0.26595 0.16229 0.185135 0.157228 0.11835 0.111086 0.13054 0.152265 0.098407 0.087037 0.119005 0.111335 0.11135 0.09954 0.091682 0.094603' and '0.0769628027961717 0.0879817673017071 0.102693260461326 0.127890558619956 0.0881286645381537 0.119938023142016 0.111385845077158 0.152741348927576 0.115265317310988 0.113149669914188 0.110216604537394 0.122942068779812 0.118612850693353 0.114615722996008 0.145983193692793 0.106001665730875 0.106819667920572 0.141533273195614 0.133994540354829 0.172064802114716 0.133244434174562 0.109426086094424 0.121738783670708 0.266476960563177 0.162579181737322 0.185465078754807 0.157354439362006 0.117599416122372 0.110415945759765 0.12958313442684 0.151131119099648 0.0979586865629627 0.0867446700508072 0.118324504684533 0.110770856198069 0.111000472772339 0.099219367161986 0.0913445613875705 0.0943421980009814' are not equal to within 0.001"
The difference between the two differs by a touch more than .001 in one case:
[1] -6.280280e-05 -1.167673e-04 -1.502605e-04 -6.255862e-05 -4.466454e-05 -1.470231e-04 -1.498451e-04 -1.163489e-04 -5.931731e-05 -6.566991e-05
[11] -5.260454e-05 -5.906878e-05 -7.785069e-05 -2.672300e-05 -5.119369e-05 -2.166573e-05 -2.066792e-05 -2.627320e-05 -1.535404e-04 -8.680211e-05
[21] -1.144342e-04 -1.308609e-05 -2.078367e-05 -5.269606e-04 -2.891817e-04 -3.300788e-04 -1.264394e-04 7.505839e-04 6.700542e-04 9.568656e-04
[31] 1.133881e-03 4.483134e-04 2.923299e-04 6.804953e-04 5.641438e-04 3.495272e-04 3.206328e-04 3.374386e-04 2.608020e-04
and may be only due to numerical precision.
Yeah, with the following modifications go ahead and check the code in: