Разреженная диагональная матрица из разреженных диагональных блоков - PullRequest
Новая
       18

Разреженная диагональная матрица из разреженных диагональных блоков

1 голос
/ 25 мая 2019

Я хочу иметь разреженную n x n матрицу A, где B - это m x m разреженная матрица и n = m * m:

Я знаю, как сделать B: 

 data = np.vstack([-np.ones(n), 4 * np.ones(n), -np.ones(n)])
 diag = np.array([-1, 0, 1])
 B = scipy.sparse.spdiags(data, diag, m, m).tocsr()

и, конечно, как сделать идентичность I.Почему-то я не могу сделать то же самое, чтобы сделать A.

Я знаю, как сделать это глупо: подставить все числа от B в A и сделать A с 5 диагоналями таким же образом, как было сделано B.Но я верю, что scipy.sparse может сделать это проще.

Я уже два часа думал об этом, может быть, найти более простой подход просто не стоит.

1 Ответ

1 голос
/ 25 мая 2019

Я обнаружил функцию sparse.block_diag. 

In [383]: n,m=3,3                                                                                        
In [385]: data = np.vstack([-np.ones(n), 4 * np.ones(n), -np.ones(n)]).astype(int) 
     ...: diag = np.array([-1, 0, 1]) 
     ...: B = sparse.spdiags(data, diag, m, m).tocsr()                                                   
In [386]: B                                                                                              
Out[386]: 
<3x3 sparse matrix of type '<class 'numpy.int64'>'
    with 7 stored elements in Compressed Sparse Row format>
In [387]: B.A                                                                                            
Out[387]: 
array([[ 4, -1,  0],
       [-1,  4, -1],
       [ 0, -1,  4]], dtype=int64)
In [388]: I = -sparse.eye(m).astype(int)

block_diag не имеет функции смещения, поэтому сначала я объединю B и I в больший блок.

In [389]: M1 = sparse.bmat([[B,I],[I,B]])                                                                
In [390]: M1                                                                                             
Out[390]: 
<6x6 sparse matrix of type '<class 'numpy.int64'>'
    with 20 stored elements in COOrdinate format>
In [391]: M1.A                                                                                           
Out[391]: 
array([[ 4, -1,  0, -1,  0,  0],
       [-1,  4, -1,  0, -1,  0],
       [ 0, -1,  4,  0,  0, -1],
       [-1,  0,  0,  4, -1,  0],
       [ 0, -1,  0, -1,  4, -1],
       [ 0,  0, -1,  0, -1,  4]], dtype=int64)
In [392]: M2 =sparse.block_diag((M1,M1))                                                                 
In [393]: M2                                                                                             
Out[393]: 
<12x12 sparse matrix of type '<class 'numpy.int64'>'
    with 40 stored elements in COOrdinate format>
In [394]: M2.A                                                                                           
Out[394]: 
array([[ 4, -1,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0],
       [-1,  4, -1,  0, -1,  0,  0,  0,  0,  0,  0,  0],
       [ 0, -1,  4,  0,  0, -1,  0,  0,  0,  0,  0,  0],
       [-1,  0,  0,  4, -1,  0,  0,  0,  0,  0,  0,  0],
       [ 0, -1,  0, -1,  4, -1,  0,  0,  0,  0,  0,  0],
       [ 0,  0, -1,  0, -1,  4,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  4, -1,  0, -1,  0,  0],
       [ 0,  0,  0,  0,  0,  0, -1,  4, -1,  0, -1,  0],
       [ 0,  0,  0,  0,  0,  0,  0, -1,  4,  0,  0, -1],
       [ 0,  0,  0,  0,  0,  0, -1,  0,  0,  4, -1,  0],
       [ 0,  0,  0,  0,  0,  0,  0, -1,  0, -1,  4, -1],
       [ 0,  0,  0,  0,  0,  0,  0,  0, -1,  0, -1,  4]], dtype=int64)

Посмотрите на код block_mat;он создает вложенный список входных данных и None и передает его в bmat.А код bmat объединяет атрибуты coo всех входов с подходящими смещениями, делая входные данные в стиле coo для матрицы результатов.sparse.hstack и sparse.vstack - другие примеры, использующие bmat.

Таким образом, следуя этим моделям, вы можете создать свой композит непосредственно из B и I.

=== *Непосредственно по диагонали 1026 *

5: 

In [405]: data = -np.ones((5,12),int)                                                                    
In [406]: data[0,:] *= -4                                                                                
In [407]: offsets=np.array([0,-1,-3,1,3])  # offsets dont have to be in order                                                         
In [408]: M4 = sparse.spdiags(data, offsets,12,12)                                                       
In [409]: M4                                                                                             
Out[409]: 
<12x12 sparse matrix of type '<class 'numpy.int64'>'
    with 52 stored elements (5 diagonals) in DIAgonal format>
In [410]: M4.A                                                                                           
Out[410]: 
array([[ 4, -1,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0],
       [-1,  4, -1,  0, -1,  0,  0,  0,  0,  0,  0,  0],
       [ 0, -1,  4, -1,  0, -1,  0,  0,  0,  0,  0,  0],
       [-1,  0, -1,  4, -1,  0, -1,  0,  0,  0,  0,  0],
       [ 0, -1,  0, -1,  4, -1,  0, -1,  0,  0,  0,  0],
       [ 0,  0, -1,  0, -1,  4, -1,  0, -1,  0,  0,  0],
       [ 0,  0,  0, -1,  0, -1,  4, -1,  0, -1,  0,  0],
       [ 0,  0,  0,  0, -1,  0, -1,  4, -1,  0, -1,  0],
       [ 0,  0,  0,  0,  0, -1,  0, -1,  4, -1,  0, -1],
       [ 0,  0,  0,  0,  0,  0, -1,  0, -1,  4, -1,  0],
       [ 0,  0,  0,  0,  0,  0,  0, -1,  0, -1,  4, -1],
       [ 0,  0,  0,  0,  0,  0,  0,  0, -1,  0, -1,  4]])

Непосредственно на coo входы: 

In [411]: data,rows,cols = [],[],[]  
In [413]: for v,z in zip([-1,-1,4,-1,-1],[-3,-1,0,1,3]): 
     ...:    len = 12-abs(z) 
     ...:    d = np.ones(len, int)*v 
     ...:    r = np.arange(len)+max(0,z) 
     ...:    c = np.arange(len)+max(0,-z) 
     ...:    data.append(d); rows.append(r); cols.append(c) 
     ...:                                                                                                
In [414]: data                                                                                           
Out[414]: 
[array([-1, -1, -1, -1, -1, -1, -1, -1, -1]),
 array([-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]),
 array([4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]),
 ...]
In [415]: rows                                                                                           
Out[415]: 
[array([0, 1, 2, 3, 4, 5, 6, 7, 8]),
 array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10]),
 ....]
In [416]: cols                                                                                           
Out[416]: 
[array([ 3,  4,  5,  6,  7,  8,  9, 10, 11]),
 array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11]),
 ...]
In [417]: M5 = sparse.coo_matrix((np.hstack(data), (np.hstack(rows), np.hstack(cols))))                  
In [418]: M5                                                                                             
Out[418]: 
<12x12 sparse matrix of type '<class 'numpy.int64'>'
    with 52 stored elements in COOrdinate format>

Или только со списком data (и смещениями): 

In [427]: M6=sparse.diags(data, [-3,-1,0,1,3],dtype=int,format='csr')                                 
In [428]: M6                                                                                          
Out[428]: 
<12x12 sparse matrix of type '<class 'numpy.int64'>'
    with 52 stored elements in Compressed Sparse Row format>
...