Manera eficiente de computar la matriz de Vandermonde.

Estoy calculando la Vandermonde matrix para una matriz 1D bastante grande. La forma natural y limpia de hacerlo es usando np.vander() . Sin embargo, encontré que esto es aprox. 2.5x más lento que un enfoque basado en la comprensión de lista.

 In [43]: x = np.arange(5000) In [44]: N = 4 In [45]: %timeit np.vander(x, N, increasing=True) 155 µs ± 205 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) # one of the listed approaches from the documentation In [46]: %timeit np.flip(np.column_stack([x**(N-1-i) for i in range(N)]), axis=1) 65.3 µs ± 235 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) In [47]: np.all(np.vander(x, N, increasing=True) == np.flip(np.column_stack([x**(N-1-i) for i in range(N)]), axis=1)) Out[47]: True 

Estoy tratando de entender dónde está el cuello de botella y la razón por la cual la implementación nativa de np.vander() es ~ 2.5x más lenta.

La eficiencia importa para mi implementación. Por lo tanto, incluso alternativas más rápidas también son bienvenidas!

Aquí hay algunos métodos más, algunos de los cuales son bastante más rápidos (en mi computadora) que los que se han publicado hasta ahora.

La observación más importante que creo es que realmente depende mucho de cuántos grados quieras. La exponencia (que creo que es especial para los pequeños exponentes enteros) solo tiene sentido para los pequeños rangos de exponentes. Cuantos más exponentes valgan mejor los enfoques basados ​​en la multiplicación.

Me gustaría resaltar un método basado en multiply.accumulate ( ma ) que es similar al enfoque integrado de numpy pero más rápido (y no porque escatimé en las comprobaciones: nnc , numpy-no-checks demuestra esto). Para todos, excepto los rangos de exponentes más pequeños, en realidad es el más rápido para mí.

Por razones que no entiendo, la implementación numpy hace tres cosas que son, a mi entender, lentas e innecesarias: (1) Hace unas cuantas copias del vector base. (2) Los hace no contiguos. (3) Hace la acumulación en el lugar que creo que obliga a amortiguar.

Otra cosa que me gustaría mencionar es que el más rápido para rangos pequeños de out_e_1 ( out_e_1 esencialmente una versión manual de ma ), se ralentiza en un factor de más de dos por la simple precaución de promocionar a un tipo de letra más grande ( safe_e_1 posiblemente un poco de un nombre inapropiado).

Los métodos basados ​​en la difusión se llaman bc_* donde * indica el eje de transmisión (b para base, e para exp) ‘truco’ significa que el resultado no es contiguo.

Tiempos (el mejor de 3):

 rep=100 n_b=5000 n_e=4 b_tp= e_tp= vander 0.16699657 ms bc_b 0.09595204 ms bc_e 0.07959786 ms ma 0.10755240 ms nnc 0.16459018 ms out_e_1 0.02037535 ms out_e_2 0.02656622 ms safe_e_1 0.04652272 ms safe_e_2 0.04081079 ms cheat bc_e_cheat 0.04668466 ms rep=100 n_b=5000 n_e=8 b_tp= e_tp= vander 0.25086462 ms bc_b apparently failed bc_e apparently failed ma 0.15843041 ms nnc 0.24713077 ms out_e_1 apparently failed out_e_2 apparently failed safe_e_1 0.15970622 ms safe_e_2 0.19672418 ms bc_e_cheat apparently failed rep=100 n_b=5000 n_e=4 b_tp= e_tp= vander 0.16225773 ms bc_b 0.53315020 ms bc_e 0.56200830 ms ma 0.07626799 ms nnc 0.16059748 ms out_e_1 0.03653416 ms out_e_2 0.04043702 ms safe_e_1 0.04060494 ms safe_e_2 0.04104209 ms cheat bc_e_cheat 0.52966076 ms rep=100 n_b=5000 n_e=8 b_tp= e_tp= vander 0.24542852 ms bc_b 2.03353578 ms bc_e 2.04281270 ms ma 0.11075758 ms nnc 0.24212880 ms out_e_1 0.14809043 ms out_e_2 0.19261359 ms safe_e_1 0.15206112 ms safe_e_2 0.19308420 ms cheat bc_e_cheat 1.99176601 ms 

Código:

 import numpy as np import types from timeit import repeat prom={np.dtype(np.int32): np.dtype(np.int64), np.dtype(float): np.dtype(float)} def RI(k, N, dt, top=100): return np.random.randint(0, top if top else N, (k, N)).astype(dt) def RA(k, N, dt, top=None): return np.add.outer(np.zeros((k,), int), np.arange(N)%(top if top else N)).astype(dt) def RU(k, N, dt, top=100): return (np.random.random((k, N))*(top if top else N)).astype(dt) def data(k, N_b, N_e, dt_b, dt_e, b_fun=RI, e_fun=RA): b = list(b_fun(k, N_b, dt_b)) e = list(e_fun(k, N_e, dt_e)) return b, e def f_vander(b, e): return np.vander(b, len(e), increasing=True) def f_bc_b(b, e): return b[:, None]**e def f_bc_e(b, e): return np.ascontiguousarray((b**e[:, None]).T) def f_ma(b, e): out = np.empty((len(b), len(e)), prom[b.dtype]) out[:, 0] = 1 np.multiply.accumulate(np.broadcast_to(b, (len(e)-1, len(b))), axis=0, out=out[:, 1:].T) return out def f_nnc(b, e): out = np.empty((len(b), len(e)), prom[b.dtype]) out[:, 0] = 1 out[:, 1:] = b[:, None] np.multiply.accumulate(out[:, 1:], out=out[:, 1:], axis=1) return out def f_out_e_1(b, e): out = np.empty((len(b), len(e)), b.dtype) out[:, 0] = 1 out[:, 1] = b out[:, 2] = c = b*b for i in range(3, len(e)): c*=b out[:, i] = c return out def f_out_e_2(b, e): out = np.empty((len(b), len(e)), b.dtype) out[:, 0] = 1 out[:, 1] = b out[:, 2] = b*b for i in range(3, len(e)): out[:, i] = out[:, i-1] * b return out def f_safe_e_1(b, e): out = np.empty((len(b), len(e)), prom[b.dtype]) out[:, 0] = 1 out[:, 1] = b out[:, 2] = c = (b*b).astype(prom[b.dtype]) for i in range(3, len(e)): c*=b out[:, i] = c return out def f_safe_e_2(b, e): out = np.empty((len(b), len(e)), prom[b.dtype]) out[:, 0] = 1 out[:, 1] = b out[:, 2] = b*b for i in range(3, len(e)): out[:, i] = out[:, i-1] * b return out def f_bc_e_cheat(b, e): return (b**e[:, None]).T for params in [(100, 5000, 4, np.int32, np.int32), (100, 5000, 8, np.int32, np.int32), (100, 5000, 4, float, np.int32), (100, 5000, 8, float, np.int32)]: k = params[0] dat = data(*params) ref = f_vander(dat[0][0], dat[1][0]) print('rep={} n_b={} n_e={} b_tp={} e_tp={}'.format(*params)) for name, func in list(globals().items()): if not name.startswith('f_') or not isinstance(func, types.FunctionType): continue try: assert np.allclose(ref, func(dat[0][0], dat[1][0])) if not func(dat[0][0], dat[1][0]).flags.c_contiguous: print('cheat', end=' ') print("{:16s}{:16.8f} ms".format(name[2:], np.min(repeat( 'f(b.pop(), e.pop())', setup='b, e = data(*p)', globals={'f':func, 'data':data, 'p':params}, number=k)) * 1000 / k)) except: print("{:16s} apparently failed".format(name[2:])) 

¿Qué hay de la exponenciación emitida?

 %timeit (x ** np.arange(N)[:, None]).T 43 µs ± 348 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

Prueba de cordura –

 np.all((x ** np.arange(N)[:, None]).T == np.vander(x, N, increasing=True)) True 

La advertencia aquí es que esta aceleración es posible solo si la matriz de entrada x tiene un dtype de int . Como @Warren Weckesser señaló en un comentario, la exponenciación transmitida se ralentiza para las matrices de punto flotante.


En cuanto a por qué np.vander es lento, eche un vistazo al código fuente :

 x = asarray(x) if x.ndim != 1: raise ValueError("x must be a one-dimensional array or sequence.") if N is None: N = len(x) v = empty((len(x), N), dtype=promote_types(x.dtype, int)) tmp = v[:, ::-1] if not increasing else v if N > 0: tmp[:, 0] = 1 if N > 1: tmp[:, 1:] = x[:, None] multiply.accumulate(tmp[:, 1:], out=tmp[:, 1:], axis=1) return v 

La función tiene que atender a muchos más casos de uso además del suyo, por lo que utiliza un método de computación más generalizado que es confiable, pero más lento (estoy apuntando específicamente a multiply.accumulate ).


Como cuestión de interés, encontré otra forma de calcular la matriz de Vandermonde, terminando con esto:

 %timeit x[:, None] ** np.arange(N) 150 µs ± 230 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

Hace lo mismo, pero es mucho más lento. La respuesta está en el hecho de que las operaciones se transmiten, pero de manera ineficiente .

Por otro lado, para los arreglos de float , esto termina siendo el mejor.