squares non marquardt levenberg leastsq least fitting python optimization numpy rotation scipy

non - levenberg marquardt python



Rotación de tensor rápido con NumPy (6)

Enfoque prospectivo y código de solución

Para la eficiencia de la memoria y, posteriormente, la eficiencia del rendimiento, podríamos usar la multiplicación de la matriz tensorial en pasos.

Para ilustrar los pasos involucrados, usemos la más simple de las soluciones con np.einsum by @pv. -

np.einsum(''ai,bj,ck,dl,abcd->ijkl'', g, g, g, g, T)

Como se ve, estamos perdiendo la primera dimensión de g con tensor-multiplicación entre sus cuatro variantes y T

Hagamos esas reducciones de suma para las multiplicaciones de la matriz tensorial en pasos. Comencemos con la primera variante de g y T :

p1 = np.einsum(''abcd, ai->bcdi'', T, g)

Por lo tanto, terminamos con un tensor de dimensiones como notación de cuerda: bcdi . Los próximos pasos implicarían reducir la suma de este tensor contra el resto de las variantes de tres g como se usa en la einsum original de einsum. Por lo tanto, la próxima reducción sería -

p2 = np.einsum(''bcdi, bj->cdij'', p1, g)

Como se ve, hemos perdido las dos primeras dimensiones con las anotaciones de cuerda: a , b . Continuamos por dos pasos más para deshacernos de d también quedará con ijkl como resultado final, como así -

p3 = np.einsum(''cdij, ck->dijk'', p2, g) p4 = np.einsum(''dijk, dl->ijkl'', p3, g)

Ahora, podríamos usar np.tensordot para estas reducciones de suma, que serían mucho más eficientes.

Implementación final

Por lo tanto, np.tensordot a np.tensordot , tendríamos la implementación final como tal -

p1 = np.tensordot(T,g,axes=((0),(0))) p2 = np.tensordot(p1,g,axes=((0),(0))) p3 = np.tensordot(p2,g,axes=((0),(0))) out = np.tensordot(p3,g,axes=((0),(0)))

Prueba de tiempo de ejecución

Probemos todos los enfoques basados ​​en NumPy publicados en otras publicaciones para resolver el problema del rendimiento.

Enfoques como funciones -

def rotT_Philipp(T, g): # @Philipp''s soln gg = np.outer(g, g) gggg = np.outer(gg, gg).reshape(4 * g.shape) axes = ((0, 2, 4, 6), (0, 1, 2, 3)) return np.tensordot(gggg, T, axes) def rotT_Sven(T, g): # @Sven Marnach''s soln Tprime = T for i in range(4): slices = [None] * 4 slices[i] = slice(None) slices *= 2 Tprime = g[slices].T * Tprime return Tprime.sum(-1).sum(-1).sum(-1).sum(-1) def rotT_pv(T, g): # @pv.''s soln return np.einsum(''ai,bj,ck,dl,abcd->ijkl'', g, g, g, g, T) def rotT_Divakar(T,g): # Posted in this post p1 = np.tensordot(T,g,axes=((0),(0))) p2 = np.tensordot(p1,g,axes=((0),(0))) p3 = np.tensordot(p2,g,axes=((0),(0))) p4 = np.tensordot(p3,g,axes=((0),(0))) return p4

Tiempos con los tamaños de conjunto de datos originales -

In [304]: # Setup inputs ...: T = np.random.rand(3,3,3,3) ...: g = np.random.rand(3,3) ...: In [305]: %timeit rotT(T, g) ...: %timeit rotT_pv(T, g) ...: %timeit rotT_Sven(T, g) ...: %timeit rotT_Philipp(T, g) ...: %timeit rotT_Divakar(T, g) ...: 100 loops, best of 3: 6.51 ms per loop 1000 loops, best of 3: 247 µs per loop 10000 loops, best of 3: 137 µs per loop 10000 loops, best of 3: 41.6 µs per loop 10000 loops, best of 3: 28.3 µs per loop In [306]: 6510.0/28.3 # Speedup with the proposed soln over original code Out[306]: 230.03533568904592

Como se discutió al comienzo de esta publicación, estamos tratando de lograr la eficiencia de la memoria y, por lo tanto, aumentar el rendimiento con ella. Vamos a probar eso a medida que aumentamos los tamaños de los conjuntos de datos:

In [307]: # Setup inputs ...: T = np.random.rand(5,5,5,5) ...: g = np.random.rand(5,5) ...: In [308]: %timeit rotT(T, g) ...: %timeit rotT_pv(T, g) ...: %timeit rotT_Sven(T, g) ...: %timeit rotT_Philipp(T, g) ...: %timeit rotT_Divakar(T, g) ...: 100 loops, best of 3: 6.54 ms per loop 100 loops, best of 3: 7.17 ms per loop 100 loops, best of 3: 2.7 ms per loop 1000 loops, best of 3: 1.47 ms per loop 10000 loops, best of 3: 39.9 µs per loop

En el corazón de una aplicación (escrita en Python y usando NumPy ) necesito rotar un tensor de 4º orden. En realidad, necesito rotar muchos tensores muchas veces y este es mi cuello de botella. Mi implementación ingenua (a continuación) que involucra ocho bucles anidados parece ser bastante lenta, pero no veo una manera de aprovechar las operaciones matriciales de NumPy y, con suerte, acelerar las cosas. Tengo la sensación de que debería estar usando np.tensordot , pero no veo cómo.

Matemáticamente, los elementos del tensor rotado, T '', están dados por: T'' ijkl = Σ g ia g jb g kc g ld T abcd con la suma sobre los índices repetidos en el lado derecho. T y Tprime son 3 * 3 * 3 * 3 matrices NumPy y la matriz de rotación g es una matriz 3 * 3 NumPy. Mi implementación lenta (tomando ~ 0.04 segundos por llamada) está debajo.

#!/usr/bin/env python import numpy as np def rotT(T, g): Tprime = np.zeros((3,3,3,3)) for i in range(3): for j in range(3): for k in range(3): for l in range(3): for ii in range(3): for jj in range(3): for kk in range(3): for ll in range(3): gg = g[ii,i]*g[jj,j]*g[kk,k]*g[ll,l] Tprime[i,j,k,l] = Tprime[i,j,k,l] + / gg*T[ii,jj,kk,ll] return Tprime if __name__ == "__main__": T = np.array([[[[ 4.66533067e+01, 5.84985000e-02, -5.37671310e-01], [ 5.84985000e-02, 1.56722231e+01, 2.32831900e-02], [ -5.37671310e-01, 2.32831900e-02, 1.33399259e+01]], [[ 4.60051700e-02, 1.54658176e+01, 2.19568200e-02], [ 1.54658176e+01, -5.18223500e-02, -1.52814920e-01], [ 2.19568200e-02, -1.52814920e-01, -2.43874100e-02]], [[ -5.35577630e-01, 1.95558600e-02, 1.31108757e+01], [ 1.95558600e-02, -1.51342210e-01, -6.67615000e-03], [ 1.31108757e+01, -6.67615000e-03, 6.90486240e-01]]], [[[ 4.60051700e-02, 1.54658176e+01, 2.19568200e-02], [ 1.54658176e+01, -5.18223500e-02, -1.52814920e-01], [ 2.19568200e-02, -1.52814920e-01, -2.43874100e-02]], [[ 1.57414726e+01, -3.86167500e-02, -1.55971950e-01], [ -3.86167500e-02, 4.65601977e+01, -3.57741000e-02], [ -1.55971950e-01, -3.57741000e-02, 1.34215636e+01]], [[ 2.58256300e-02, -1.49072770e-01, -7.38843000e-03], [ -1.49072770e-01, -3.63410500e-02, 1.32039847e+01], [ -7.38843000e-03, 1.32039847e+01, 1.38172700e-02]]], [[[ -5.35577630e-01, 1.95558600e-02, 1.31108757e+01], [ 1.95558600e-02, -1.51342210e-01, -6.67615000e-03], [ 1.31108757e+01, -6.67615000e-03, 6.90486240e-01]], [[ 2.58256300e-02, -1.49072770e-01, -7.38843000e-03], [ -1.49072770e-01, -3.63410500e-02, 1.32039847e+01], [ -7.38843000e-03, 1.32039847e+01, 1.38172700e-02]], [[ 1.33639532e+01, -1.26331100e-02, 6.84650400e-01], [ -1.26331100e-02, 1.34222177e+01, 1.67851800e-02], [ 6.84650400e-01, 1.67851800e-02, 4.89151396e+01]]]]) g = np.array([[ 0.79389393, 0.54184237, 0.27593346], [-0.59925749, 0.62028664, 0.50609776], [ 0.10306737, -0.56714313, 0.8171449 ]]) for i in range(100): Tprime = rotT(T,g)

¿Hay alguna manera de hacer que esto vaya más rápido? Hacer que el código se generalice a otros rangos de tensor sería útil, pero es menos importante.


Aquí es cómo hacerlo con un solo bucle de Python:

def rotT(T, g): Tprime = T for i in range(4): slices = [None] * 4 slices[i] = slice(None) slices *= 2 Tprime = g[slices].T * Tprime return Tprime.sum(-1).sum(-1).sum(-1).sum(-1)

Es cierto que esto es un poco difícil de entender a primera vista, pero es bastante más rápido :)


Gracias al arduo trabajo de M. Wiebe, la próxima versión de Numpy (que probablemente será 1.6) lo hará aún más fácil:

>>> Trot = np.einsum(''ai,bj,ck,dl,abcd->ijkl'', g, g, g, g, T)

Sin embargo, el enfoque de Philipp es en este momento 3 veces más rápido, pero tal vez haya algún margen de mejora. La diferencia de velocidad probablemente se deba principalmente a que el tensordot puede desenrollar toda la operación como un solo producto de matriz que puede pasarse a BLAS, y así evitar la sobrecarga asociada con arreglos pequeños --- esto no es posible para Einstein general suma, ya que no todas las operaciones que se pueden expresar en esta forma se resuelven en un único producto de matriz.


Para usar el tensordot , calcule el producto externo de los tensores g :

def rotT(T, g): gg = np.outer(g, g) gggg = np.outer(gg, gg).reshape(4 * g.shape) axes = ((0, 2, 4, 6), (0, 1, 2, 3)) return np.tensordot(gggg, T, axes)

En mi sistema, esto es alrededor de siete veces más rápido que la solución de Sven. Si el tensor g no cambia con frecuencia, también puede almacenar en caché el tensor gggg . Si hace esto y tensordot ( tensordot código de tensordot , sin controles, sin formas genéricas), puede hacerlo dos veces más rápido:

def rotT(T, gggg): return np.dot(gggg.transpose((1, 3, 5, 7, 0, 2, 4, 6)).reshape((81, 81)), T.reshape(81, 1)).reshape((3, 3, 3, 3))

Resultados del timeit en mi computadora portátil hogareña (500 iteraciones):

Your original code: 19.471129179 Sven''s code: 0.718412876129 My first code: 0.118047952652 My second code: 0.0690279006958

Los números en mi máquina de trabajo son:

Your original code: 9.77922987938 Sven''s code: 0.137110948563 My first code: 0.0569641590118 My second code: 0.0308079719543


Pensé que contribuiría con un punto de datos relativamente nuevo a estos puntos de referencia, utilizando parakeet , uno de los compiladores de JIT nudosos que surgieron en los últimos meses. (El otro que conozco es numba , pero no lo numba aquí).

Después de superar el proceso de instalación un tanto laberíntico para LLVM, puede decorar muchas funciones numpy para (a menudo) acelerar su rendimiento:

import numpy as np import parakeet @parakeet.jit def rotT(T, g): # ...

Solo intenté aplicar el JIT al código de Andrew en la pregunta original, pero funciona bastante bien (> 10 veces la aceleración) por no tener que escribir ningún código nuevo:

andrew 10 loops, best of 3: 206 msec per loop andrew_jit 10 loops, best of 3: 13.3 msec per loop sven 100 loops, best of 3: 2.39 msec per loop philipp 1000 loops, best of 3: 0.879 msec per loop

Para estos tiempos (en mi computadora portátil) ejecuté cada función diez veces, para dar al JIT la oportunidad de identificar y optimizar las rutas de los códigos de acceso rápido.

Curiosamente, ¡las sugerencias de Sven y Philipp siguen siendo órdenes de magnitud más rápidas!


Por curiosidad, Cython implementación de Cython de un código ingenuo de la pregunta con el código numpy de la respuesta de @ Philipp . El código de Cython es cuatro veces más rápido en mi máquina:

#cython: boundscheck=False, wraparound=False import numpy as np cimport numpy as np def rotT(np.ndarray[np.float64_t, ndim=4] T, np.ndarray[np.float64_t, ndim=2] g): cdef np.ndarray[np.float64_t, ndim=4] Tprime cdef Py_ssize_t i, j, k, l, ii, jj, kk, ll cdef np.float64_t gg Tprime = np.zeros((3,3,3,3), dtype=T.dtype) for i in range(3): for j in range(3): for k in range(3): for l in range(3): for ii in range(3): for jj in range(3): for kk in range(3): for ll in range(3): gg = g[ii,i]*g[jj,j]*g[kk,k]*g[ll,l] Tprime[i,j,k,l] = Tprime[i,j,k,l] + / gg*T[ii,jj,kk,ll] return Tprime