Page 1 of 1

Cholesky decomposition question

PostPosted: Mon Nov 15, 2010 3:03 pm
by levgivon
I'm observing errors in the output of culaDeviceDpotrf; specifically, it seems that the off-diagonal values in the decomposition are not problematic. Decomposing the following positive definite matrix, for instance,

Code: Select all
array([[ 0.03691047,  0.12910772,  0.08621521,  0.16356607],
           [ 0.12910772,  0.99419601,  0.36184852,  1.24055282],
           [ 0.08621521,  0.36184852,  0.40656805,  0.49676872],
           [ 0.16356607,  1.24055282,  0.49676872,  1.56339406]])

results in the upper diagonal matrix

Code: Select all
array([[ 0.19212098,  0.12910772,  0.08621521,  0.16356607],
           [ 0.                ,  0.73661053,  0.36184852,  1.24055282],
           [ 0.                ,  0.                 ,  0.44552265,  0.49676872],
           [ 0.                ,  0.                 ,  0.                ,  0.08301626]])

In contrast, performing the decomposition on the CPU with Scientific Python results in a correct answer:

Code: Select all
array([[ 0.19212098,  0.6720126  ,  0.44875477,  0.85137012],
           [ 0.                ,  0.73661053,  0.08183384,  0.90742847],
           [ 0.                ,  0.                ,  0.44552265,  0.09080114],
           [ 0.                ,  0.                ,  0.                ,  0.08301626]])


What could be causing the error?
I've attached a Python program that demonstrates the problem; it uses pycuda and scipy. Also, I'm using the premium release of CULA 2.1 on 64-bit Linux.

Re: Cholesky decomposition question

PostPosted: Mon Nov 15, 2010 4:43 pm
by kyle
I just quickly ran your test matrix through CULA (using a simple MATLAB script) and got the answer that matches your "correct answer".

I'm not Python or PyCUDA expert, but I'd suggest double checking that code to see if something is wrong there. Alternatively, check out the PyCULA work that has been done if you'd like to easily call CULA from a Python environment.

Code: Select all
>> R = culachol(A)

R =

    0.1921    0.6720    0.4488    0.8514
         0    0.7366    0.0818    0.9074
         0         0    0.4455    0.0908
         0         0         0    0.0830

>> norm(A - R*R')

ans =

1.1102e-016

Re: Cholesky decomposition question

PostPosted: Mon Nov 15, 2010 6:22 pm
by levgivon
Figured out the problem; I wasn't taking into consideration CULA's assumptions regarding matrix storage.