Estoy buscando una forma elegante de calcular el "producto" de una convolución discreta, en lugar de la suma.

Aquí está la fórmula de una convolución discreta :

enter image description here

En este caso podemos usar: conv(x,y)

Ahora me gustaría implementar esas operaciones

enter image description here

Por supuesto que puedo usar un bucle, pero estoy buscando un truco para linealizar esta operación.

EJEMPLO:

f = [2 4 3 9 7 1]
g = [3 2 1]

dist = length(g)-1;

for ii = 1:length(f)-dist
    x(ii) = prod(f(ii:ii+dist).*g)
end

x =

144    648   1134    378
4
obchardon 28 dic. 2016 a las 20:16

3 respuestas

La mejor respuesta

Otra solución inspirada en parte por respuesta Dev-iL a relativamente la misma pregunta

exp(sum(log(g))+conv(log(f),[1 1 1],'valid'))

O

exp(sum(log(g))+movsum(log(f),numel(g),'Endpoints', 'discard'))

Desde exp(sum(log(x))) = prod(x)

Pero aquí, en lugar de un vector, tenemos dos vectores f y g.

La fórmula deseada se puede reformular como:

enter image description here

Tiempo en octava:

f= rand(1,1000000);
g= rand(1,100);

disp('----------EXP(LOG)------')
tic
    exp(sum(log(g))+conv(log(f),ones(1,numel(g))));
toc

disp('----------BSXFUN------')
tic
    ind = bsxfun(@plus, 0:numel(f)-numel(g), (1:numel(g)).');
    x = prod(bsxfun(@times, f(ind), g(:)),1);
toc
disp('----------HANKEL------')
tic
    prod(g)*prod(hankel(f(1:numel(g)), f(numel(g):end)));
toc

disp('----------CUMPROD-----')
tic
    pf = cumprod(f);
    x = prod(g)*pf(numel(g):end)./[1 pf(1:(end-numel(g)))];
toc

Resultado:

----------EXP(LOG)------%rahnema1
Elapsed time is 0.211445 seconds.
----------BSXFUN--------%Luis Mendo
Elapsed time is 1.94182 seconds.
----------HANKEL--------%gnovice
Elapsed time is 1.46593 seconds.
----------CUMPROD-------%gnovice
Elapsed time is 0.00748992 seconds.
4
Community 23 may. 2017 a las 12:16

Aquí hay una forma de evitar bucles:

ind = bsxfun(@plus, 0:numel(f)-numel(g), (1:numel(g)).');
x = prod(bsxfun(@times, f(ind), g(:)), 1);

Esto funciona de la siguiente manera. ind representa la ventana deslizante de índices, cada columna corresponde a un desplazamiento. Por ejemplo, si g tiene el tamaño 3, la matriz ind será

1 2 3 4 ...
2 3 4 5 ...
3 4 5 6 ...

Esto se usa para indexar en f. El resultado se multiplica (con difusión) por g como una columna. Finalmente se calcula el producto de los elementos de cada columna.

2
Luis Mendo 28 dic. 2016 a las 18:16

Solución cumprod: (muy eficiente)

>> pf = cumprod(f);
>> x = prod(g).*pf(numel(g):end)./[1 pf(1:(end-numel(g)))]

x =

         144         648        1134         378

Este primero toma el producto acumulativo de f usando cumprod. Al dividir cada elemento por el producto acumulativo 3 elementos anteriores, obtenemos el producto de cada numel(g) - amplia ventana deslizante a lo largo de f. Luego simplemente multiplique por el producto de los elementos de g.

NOTA: cuando f tiene muchos elementos, o valores extremos (grandes o pequeños), puede encontrarse con problemas con la precisión o con el desbordamiento / desbordamiento al realizar el producto acumulativo. Una posible forma de mitigar esto sería aplicar una escala a f antes del producto acumulativo y luego deshacerlo:

c = ...set a scaling factor...
pf = cumprod(f./c);
x = prod(c.*g).*pf(numel(g):end)./[1 pf(1:(end-numel(g)))];

La opción para c podría ser algo como mean(abs(f)) o max(abs(f)) para que el f escalado dé resultados que estén mejor delimitados (es decir, valores más cercanos a 1). Esto no cambia apreciablemente los resultados de tiempo a continuación.


Solución hankel: (no tan eficiente pero aún interesante)

>> x = prod(g).*prod(hankel(f(1:numel(g)), f(numel(g):end)))

x =

         144         648        1134         378

La llamada a hankel crea una matriz donde cada columna tiene el contenido de una de las numel(g) - ventanas deslizantes anchas. Tomar el producto en cada columna y luego multiplicarlo por el producto de los elementos de g da su respuesta. Sin embargo, para los vectores grandes f y / o g esto podría implicar una gran cantidad de cómputo adicional y utilizar mucha memoria.


Momento de los resultados:

Probé 6 soluciones (el bucle en su pregunta, 2 soluciones de rahnema (conv(log) y {{X1} }), la solución bsxfun de Luis, y mi cumprod y hankel soluciones) usando f = rand(1,1000000); y g = rand(1,100); y promediando más de 40 iteraciones. Esto es lo que obtuve (ejecutando Windows 7 x64, 16 GB de RAM, MATLAB R2016b):

solution    | avg. time (s)
------------+---------------
loop        |       1.10671
conv(log)   |       0.04984
movsum(log) |       0.03736
bsxfun      |       1.20472
cumprod     |       0.01469
hankel      |       1.17704
5
Community 23 may. 2017 a las 12:24