7

As a small exercise before i start playing with numeric code in python I am trying to make an LDLT algorithm. Just to "get the feet wet".

However I seem to be lacking a fundamental understanding of the numpy array. See the following example:

def ldlt(Matrix):
    import numpy

    (NRow, NCol) = Matrix.shape

    for col in range(NCol):
        Tmp = 1/Matrix[col,col]
        for D in range(col+1, NCol):
            Matrix[col,D] = Matrix[D,col]*Tmp  
            
if __name__ == '__main__':
    import numpy
    A = numpy.array([[2,-1,0],[-1,2,-1],[0,-1,2]])
    ldlt(A)

The example is not the full code I am working on. However, try and run it, and set a break-point at Matrix[col,D] = ...

What I expect for the first evaluation is that row 0 column 1 (starting value of -1) to be set equal to = -1*(1/2) = -0.5.

However when running the code it seems to be set equal to 0. Why ? There must be something fundamental which I have not really understood?

Thanks in advance for all of you guys helping me out.

EDIT 1:

Python Ver.: 3.3 Tmp.: become 0.5 (As reported by my debugger).

2 Answers 2

4

The following may show what's going on:

>>> A = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])
>>> A.dtype
dtype('int32')
>>> A[0, 1]
-1
>>> A[0, 1] * 0.5
-0.5
>>> A[0, 1] *= 0.5
>>> A[0, 1]
0
>>> int(-0.5)
0

Your array can only hold 32-bit integers, so any floating point value you try to assign to it will be cast, i.e. truncated, to an int32.


For the same price, here's a more numpythonic way of doing what you were after: for loops are generally to be avoided, as they defeat the whole purpose of numpy:

def ldlt_np(arr) :
    rows, cols = arr.shape
    tmp = 1 / np.diag(arr) # this is a float array
    mask = np.tril_indices(cols)
    ret = arr * tmp[:, None] # this will also be a float array
    ret[mask] = arr[mask]

    return ret

>>> A = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])
>>> ldlt_np(A)
array([[ 2. , -0.5,  0. ],
       [-1. ,  2. , -0.5],
       [ 0. , -1. ,  2. ]])
Sign up to request clarification or add additional context in comments.

1 Comment

Yup that was it. Thanks a lot - It would have been a long time before i had stumbled across that.
0

numpy arrays have fixed type. You can't change an int array to floats later. Initialize the array as an array of floats:

A = numpy.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]], numpy.float)

1 Comment

I am on python 3.3 . Checking my piece of code - Tmp become 0.5 at evaluation.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.