Cholesky Matrix Function

4 replies [Last post]
neale's picture
Offline
Joined: 07/31/2009

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);
}
}

neale's picture
Offline
Joined: 07/31/2009
Some progress but I'm missing a step - help?

I did the following:

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...

mspiegel's picture
Offline
Joined: 07/31/2009
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:

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 backend

neale's picture
Offline
Joined: 07/31/2009
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 :
[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.

mspiegel's picture
Offline
Joined: 07/31/2009
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.