Cholesky Matrix Function
neale
Joined: 07/31/2009
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);
}
}
Some progress but I'm missing a step - help?
1. Added to MxAlgebraFunctions.R
cholesky <- function(x) {
x <- as.matrix(x)
if(nrow(x) != ncol(x)) {
stop("matrix must be square")
}
retval <- .Call("omxCallAlgebra",
list(x), # Flatten args into a single list
imxLookupSymbolTable("cholesky"),
NA)
return(as.matrix(as.numeric(retval)))
}
2. Added to the end of omxSymbolTable.tab
57 "Cholesky Decomposition" "cholesky" 1 "omxCholesky"
3. Added to omxAlgebraFunctions.c
void omxCholesky(omxMatrix** matList, int numArgs, omxMatrix* result)
{
if(OMX_DEBUG_ALGEBRA) { Rprintf("ALGEBRA: Matrix Cholesky Decomposition.\n");}
omxMatrix* inMat = matList[0];
int l = 0; char u = 'U';
omxCopyMatrix(result, inMat);
if(result->rows != result->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);
}
F77_CALL(dpotrf)(&u, &(result->rows), result->data, &(result->cols), &l);
if(l != 0) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Attempted to Cholesky decompose non-invertable matrix.\n");
omxRaiseError(result->currentState, -1, errstr);
Free(errstr);
}
}
(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...
Log in or register to post comments
In reply to Some progress but I'm missing a step - help? by neale
On second thought, is it true
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:2. Added to the end of omxSymbolTable.tab
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:
A <- mxMatrix('Full', 4, 4, values = foo, 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 backendLog in or register to post comments
In reply to On second thought, is it true by mspiegel
Ok got it
>
> 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.
Log in or register to post comments
In reply to Ok got it by neale
Yeah, with the following
Log in or register to post comments