List as index in matlab

Question

There is this method written in Matlab that I want to translate into Python. However, I don't understand how to interpret the notation of indexing the sparse matrix M with a row of the matrix faces. What would be the equivalent in Python?

M = spalloc(size(template,1), size(template,1), 10*size(template,1));

for i = 1:size(faces,1)
    v = faces(i,:); % faces is a Nx3 matrix
    ...

    M(v, v) = M(v, v) + WIJ; % WIJ is some 3x3 matrix

Answer 1

@Eric Yu` uses a dense numpy array:

In [239]: A=np.array([[1,2,3],[3,4,5],[5,6,7]])                              
In [240]: A                                                                  
Out[240]: 
array([[1, 2, 3],
       [3, 4, 5],
       [5, 6, 7]])
In [241]: v=[0,1]                                                            

this indexing selects rows:

In [242]: A[v]                                                               
Out[242]: 
array([[1, 2, 3],
       [3, 4, 5]])

and from that select columns:

In [243]: A[v][:,v]                                                          
Out[243]: 
array([[1, 2],
       [3, 4]])

But A[v] is a copy, not a view, so assignment will fail:

In [244]: A[v][:,v] = 0                                                      
In [245]: A                                                                  
Out[245]: 
array([[1, 2, 3],
       [3, 4, 5],
       [5, 6, 7]])

===

To properly index a block of a numpy array, use ix_ (or equivalent) to create indexing arrays that broadcast against each other to define the block:

In [247]: np.ix_(v,v)                                                        
Out[247]: 
(array([[0],
        [1]]), array([[0, 1]]))
In [248]: A[np.ix_(v,v)]                                                     
Out[248]: 
array([[1, 2],
       [3, 4]])
In [249]: A[np.ix_(v,v)]=0                                                   
In [250]: A                                                                  
Out[250]: 
array([[0, 0, 3],
       [0, 0, 5],
       [5, 6, 7]])

Without the ix_ transform, indexing with [v,v] selects a diagonal:

In [251]: A[v,v]                                                             
Out[251]: array([0, 0])

MATLAB M(v,v) indexes the block. Indexing the diagonal on the other hand requires use of sub2idx (or something like that). This is a case where MATLAB's indexing notation makes one task easy, and the other more complex. numpy does the reverse.

===

What I wrote is applicable to sparse matrices as well

In [253]: M=sparse.lil_matrix(np.array([[1,2,3],[3,4,5],[5,6,7]]))           
In [254]: M                                                                  
Out[254]: 
<3x3 sparse matrix of type '<class 'numpy.int64'>'
    with 9 stored elements in LInked List format>

The diagonal selection:

In [255]: M[v,v]                                                             
Out[255]: 
<1x2 sparse matrix of type '<class 'numpy.int64'>'
    with 2 stored elements in LInked List format>
In [256]: _.A                                                                
Out[256]: array([[1, 4]], dtype=int64)

Note that this matrix is (1,2), still 2d, in the style of MATLAB matrices.

block selection:

In [258]: M[np.ix_(v,v)]                                                     
Out[258]: 
<2x2 sparse matrix of type '<class 'numpy.int64'>'
    with 4 stored elements in LInked List format>
In [259]: _.A                                                                
Out[259]: 
array([[1, 2],
       [3, 4]], dtype=int64)
In [260]: M[np.ix_(v,v)]=0                                                   
In [261]: M.A                                                                
Out[261]: 
array([[0, 0, 3],
       [0, 0, 5],
       [5, 6, 7]], dtype=int64)

sparse.csr_matrix will index in the same way (with some differences in the assignment step).

Answer 2

import numpy as np
A=[[1,2,3],[3,4,5],[5,6,7]]
M=np.array(A)
v=[0,1]
M[v][:,v]

the result is:

array([[1, 2],
       [3, 4]])