Search code examples
pythonarraysnumpyvariable-assignmentoperation

Basic Numpy array value assignment


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


Solution

  • 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. ]])