Dove spfun incontra bsxfun

Question

Esiste un modo per eseguire operazioni element-like come bsxfun, ma solo funziona sugli elementi diversi da zero di una matrice sparsa?

In particolare, per ciascun elemento diverso da zero nella matrice A alla posizione (i,j) voglio trovare il prodotto di tutti gli elementi diversi da zero nella fila -esima i, tranne elemento (i,j).

Per esempio se la riga i-th assomiglia a questo:

0 5 3 0 0 4 0 0 

Il risultato dovrebbe essere simile a questo:

0 12 20 0 0 15 0 0 

La soluzione più ovvia sembra essere quello di prendere il prodotto della non zero elementi lungo ciascuna riga e quindi dividere ciascun elemento indietro rispetto al prodotto della riga. Quindi, nell'esempio sopra, il prodotto di riga è 5 x 3 x 4 = 60 e quindi ho appena diviso 53 e 4 nelle rispettive posizioni.

Data una matrice sparsa A, questa è la mia soluzione migliore finora:

[M N] = size(A); 
[row col vals] = find(A); 
row_prod = accumarray(row,vals,[],@prod); 
B = bsxfun(@ldivide, A, row_prod); 
B = sparse(row,col,B(sub2ind([M N],row,col)),M,N); 

Le prime tre righe realizzare ciò che voglio: un vettore colonna rappresenta il prodotto degli elementi non nulli di ogni riga. Tuttavia, ci sono alcuni problemi con le ultime due righe.

bsxfun sta per restituire una matrice non-sparse la dimensione di A

Si sta per perdere un sacco di cicli di divisione per zero inutilmente.

Il risultato è una matrice composta principalmente da Inf o -Inf, in cui desidero solo zeri.

Non riesco a mascherare lo Inf s perché Matlab definisce zero volte l'infinito per essere NaN.

Ho solo bisogno di mordere il proiettile e scrivere un ciclo for per questo? O c'è un altro modo per affrontarlo?

Answer 1

Penso di aver trovato una soluzione che risolve la maggior parte delle mie preoccupazioni di cui sopra. A volte quando ho il martello bsxfun nella mia mano, il mondo intero inizia a sembrare un chiodo. Dimentico che alcune delle semplici moltiplicazioni che faccio con bsxfun potrebbero essere risolte altrettanto facilmente (e probabilmente più facilmente) usando la moltiplicazione di matrici. Anche se non penso che la mia soluzione a questo problema sia più leggibile, è un ordine di grandezza più efficiente della mia ultima soluzione.

% 'This happens once, outside the loop, since the size' 
% 'and sparsity structure of A dont change in the loop' 
[M N] = size(A); 
[row col] = find(A); 

%% 'Inside iterative loop' 

    % 'Get the product of non-zero row elements' 
    row_prod = accumarray(row,nonzeros(A),[],@prod); 

    % 'Use row products to compute 'leave-one-out' row products' 
    B = spdiags(row_prod,0,M,M)*spfun(@(x) 1./x, A); 

Sarei ancora interessato a sentire altri suggerimenti se questo può essere migliorato.