J'utilise scipy pour générer une matrice de différence finie clairsemée, la construisant initialement à partir de matrices de blocs, puis éditant la diagonale pour tenir compte des conditions aux limites. La matrice clairsemée résultante est du type BSR. J'ai trouvé que si je convertis la matrice en une matrice dense puis en une matrice clairsemée en utilisant la fonction scipy.sparse.BSR_matrix, je me retrouve avec une matrice plus clairsemée qu'auparavant. Voici le code que j'utilise pour générer la matrice:

size = (4,4)

xDiff = np.zeros((size[0]+1,size[0]))
ix,jx = np.indices(xDiff.shape)
xDiff[ix==jx] = 1
xDiff[ix==jx+1] = -1

yDiff = np.zeros((size[1]+1,size[1]))
iy,jy = np.indices(yDiff.shape)
yDiff[iy==jy] = 1
yDiff[iy==jy+1] = -1

Ax = sp.sparse.dia_matrix(-np.matmul(np.transpose(xDiff),xDiff))
Ay = sp.sparse.dia_matrix(-np.matmul(np.transpose(yDiff),yDiff))

lap = sp.sparse.kron(sp.sparse.eye(size[1]),Ax) + sp.sparse.kron(Ay,sp.sparse.eye(size[0]))

#set up boundary conditions
BC_diag = np.array([2]+[1]*(size[0]-2)+[2]+([1]+[0]*(size[0]-2)+[1])*(size[1]-2)+[2]+[1]*(size[0]-2)+[2])

lap += sp.sparse.diags(BC_diag)

Si je vérifie la rareté de cette matrice, je vois ce qui suit:

lap
<16x16 sparse matrix of type '<class 'numpy.float64'>'
with 160 stored elements (blocksize = 4x4) in Block Sparse Row format>

Cependant, si je le convertis en une matrice dense, puis reviens au même format clairsemé, je vois une matrice beaucoup plus clairsemée:

sp.sparse.bsr_matrix(lap.todense())
<16x16 sparse matrix of type '<class 'numpy.float64'>'
with 64 stored elements (blocksize = 1x1) in Block Sparse Row format>

Je soupçonne que cela se produit parce que j'ai construit la matrice à l'aide de la fonction sparse.kron, mais ma question est de savoir s'il existe un moyen d'arriver à la plus petite matrice clairsemée sans convertir d'abord en dense, par exemple si je finis par vouloir simuler un très grand domaine.

3
algol 14 mars 2019 à 21:04

2 réponses

Meilleure réponse

BSR stocke les données dans des blocs denses:

In [167]: lap.data.shape                                                        
Out[167]: (10, 4, 4)

Dans ce cas, ces blocs ont plusieurs zéros.

In [168]: lap1 = lap.tocsr() 
In [170]: lap1                                                                  
Out[170]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 160 stored elements in Compressed Sparse Row format>
In [171]: lap1.data                                                             
Out[171]: 
array([-2.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,  1., -3.,  1.,  0.,  0.,
        1.,  0.,  0.,  0.,  1., -3.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,
        1., -2.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0., -3.,  1.,  0.,
        0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  1., -4.,  1.,  0., 
        ...
        0.,  0.,  1., -2.])

Nettoyage en place:

In [172]: lap1.eliminate_zeros()                                                
In [173]: lap1                                                                  
Out[173]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 64 stored elements in Compressed Sparse Row format>

Si je spécifie le format csr lors de l'utilisation de kron:

In [181]: lap2 = sparse.kron(np.eye(size[1]),Ax,format='csr') + sparse.kron(Ay,n
     ...: p.eye(size[0]), format='csr')                                         
In [182]: lap2                                                                  
Out[182]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 64 stored elements in Compressed Sparse Row format>
2
hpaulj 14 mars 2019 à 19:02

[J'ai été informé que ma réponse est incorrecte. La raison, si je comprends bien, est que Scipy n'utilise pas Lapack pour créer des matrices mais utilise son propre code à cet effet. Intéressant. L'information, bien qu'inattendue, a le pouvoir. Je m'en remettrai!

[Je laisserai la réponse affichée pour référence, mais je n'affirmerai plus que la réponse était correcte.]

De manière générale, en ce qui concerne les structures de données complexes comme les matrices clairsemées, vous avez deux cas:

  1. le constructeur connaît à l'avance le contenu complet de la structure; ou
  2. la structure est conçue pour être construite progressivement afin que le contenu complet de la structure ne soit connu qu'une fois la structure terminée.

Le cas classique de la structure de données compliquée est le cas de l'arbre binaire. Vous pouvez rendre un arbre binaire plus efficace en le copiant une fois qu'il est terminé. Sinon, l'implémentation standard rouge-noir de l'arborescence laisse certains chemins de recherche deux fois plus longs que d'autres, ce qui est généralement correct mais n'est pas optimal.

Maintenant, vous saviez probablement tout cela, mais je le mentionne pour une raison. Scipy dépend de Lapack. Lapack apporte plusieurs schémas de stockage différents. Deux d'entre eux sont les

  • clairsemée et
  • bagué

Régimes. Il semblerait que Scipy commence par stocker votre matrice comme clairsemée, où les indices de chaque élément non nul sont explicitement stockés; mais que, sur copie, Scipy remarque que la représentation en bandes est la plus appropriée - car votre matrice est, après tout, en bandes.

-1
thb 14 mars 2019 à 19:41