# Cholesky Matrix Function

5 posts / 0 new Offline
Joined: 07/31/2009 - 15:14
Cholesky Matrix Function

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;

int ipiv[inMat-&gt;rows], lwork;
lwork = 4 * inMat-&gt;rows * inMat-&gt;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-&gt;currentState, -1, errstr);
Free(errstr);
}

omxCopyMatrix(result, inMat);
F77_CALL(dpbtrf)('L', &amp;(result-&gt;cols), result-&gt;data, &amp;(result-&gt;leading), ipiv, &amp;l);
if(l != 0) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Attempted to Choleskify non-invertable matrix.");
omxRaiseError(result-&gt;currentState, -1, errstr);
Free(errstr);
}


} Offline
Joined: 07/31/2009 - 15:14
Some progress but I'm missing a step - help?

I did the following:

cholesky <- function(x) {
x <- as.matrix(x)
if(nrow(x) != ncol(x)) {
stop("matrix must be square")
}

retval &lt;- .Call("omxCallAlgebra",
list(x),         # Flatten args into a single list
imxLookupSymbolTable("cholesky"),
NA)

return(as.matrix(as.numeric(retval)))


}

1. Added to the end of omxSymbolTable.tab

57 "Cholesky Decomposition" "cholesky" 1 "omxCholesky"

void omxCholesky(omxMatrix** matList, int numArgs, omxMatrix* result)
{

if(OMX_DEBUG_ALGEBRA) { Rprintf("ALGEBRA: Matrix Cholesky Decomposition.\n");}

omxMatrix* inMat = matList;

int l = 0; char u = 'U';
omxCopyMatrix(result, inMat);
if(result-&gt;rows != result-&gt;cols) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Cholesky decomposition of non-square matrix cannot even be attempted\n");
omxRaiseError(result-&gt;currentState, -1, errstr);
Free(errstr);
}

F77_CALL(dpotrf)(&amp;u, &amp;(result-&gt;rows), result-&gt;data, &amp;(result-&gt;cols), &amp;l);
if(l != 0) {
char *errstr = Calloc(250, char);
sprintf(errstr, "Attempted to Cholesky decompose non-invertable matrix.\n");
omxRaiseError(result-&gt;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... Offline
Joined: 07/31/2009 - 15:24
On second thought, is it true

On second thought, is it true that you want the behavior of cholesky(A) to be identical to calling the R function chol(A)?

If that's the case, then don't define the function cholesky(A) in R. Just use the R function chol(A), and then define the function omxCholesky in C so that the back-end knows what to do. I would change your directions to:

1. 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 backend Offline
Joined: 07/31/2009 - 15:14
Ok got it

> 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 :
 "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:

 -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
 -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
 -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
 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. Offline
Joined: 07/31/2009 - 15:24
Yeah, with the following

Yeah, with the following modifications go ahead and check the code in:

• Add the test case(s) to models/passing/AlgebraComputePassing.R. That file contains the battery of tests on our matrix algebra routines.
• Increase the error tolerance on the failing test to 0.05 or whatever is necessary. The IntroSEM-Foo.R tests sometimes wiggle around.