Search code examples
pythonmatlabmatrixsparse-matrix

Python correspondent for MATLAB matrix operation


I have a vector of indices (let's call it peo), a sparse matrix P and a matrix W. In MATLAB I can do an operation of this kind:

P(peo, peo) = W(peo, peo)

Is there a way to do the same in Python maintaining the same computational and time complexity?


Solution

  • Choosing library

    There is a very similar way of expressing that in python if you use dense matrices. Using a sparse matrix is a little more complex. In general, if your code is not slowed by dense matrices too much and memory is not a problem I would stick to dense matrices with numpy since it is very convenient. (As they say premature optimization is the root of all evil... or something like that). However if you really need sparse matrices scipy will offer you an option for that.

    Dense matrices

    If you want to use dense matrices, you can use numpy to define matrices and peo should be defined as a list. Then you should note that the indexing with list (or vectors) doesn't work the same way in matlab as it does in python. Check this for details. (thank you Cris Lunego for pointing that out). To circumvent this and obtain the same behaviour as matlab, we will be using numpy.ix_. This will allows us to reproduce the matlab behavior of indexing with minimal alteration to our code.

    Here is an example:

    import numpy as np
    
    # Dummy matrices definition
    peo = [1, 3, 4]
    P = np.zeros((5, 5))
    W = np.ones((5, 5))
    
    # Assignment (using np.ix_ to reproduce matlab behavior)
    P[np.ix_(peo, peo)] = W[np.ix_(peo, peo)]  
    
    print(P)
    

    Sparse matrices

    For sparse matrices, scipy has a package called sparse that allows you to use sparse matrices along the way of matlab does. It gives you an actual choice on how the matrix should be represented where matlab doesn't. With great power comes great responsabiliy. Taking the time to read the pros and cons of each representation will help you choose the right one for your application.

    In general it's hard to guarantee the exact same complexity because the languages are different and I don't know the intricate details of each. But the concept of sparse matrices is the same in scipy and matlab so you can expect the complexity to be comparable. (You might even be faster in python since you can choose a representation tailored to your needs).

    Note that in this case if you want to keep working the same way as you describe in matlab you should choose a dok or lil representation. Those are the only two formats that allow efficient index access and sparsity change.

    Here is an example of what you want to archive using dok representation:

    from scipy.sparse import dok_matrix
    import numpy as np
    
    
    # Dummy matrices definition
    peo = [1, 2, 4]
    P = dok_matrix((5, 5))
    W = np.ones((5, 5))
    
    # Assignment (using np.ix_ to reproduce matlab behavior)
    P[np.ix_(peo, peo)] = W[np.ix_(peo, peo)]
    
    print(P.toarray())
    

    If you are interested in the pros and cons of sparse matrix representation and algebra in Python here is a post that explores this a bit as well as performances. This is to take with a grain of salt since it is a little old, but the ideas behind it are still mostly correct.