python - transpuesta - numpy tutorial español pdf

Multiplicación de matrices de enteros en NumPy/SciPy

Question

Haciendo algo como

import numpy as np
a = np.random.rand(10**4, 10**4)
b = np.dot(a, a)

utiliza múltiples núcleos, y se ejecuta muy bien.

Sin embargo, los elementos en a son flotadores de 64 bits (o 32 bits en plataformas de 32 bits?), Y me gustaría multiplicar las matrices de enteros de 8 bits. Intentando lo siguiente, sin embargo:

a = np.random.randint(2, size=(n, n)).astype(np.int8)

hace que el producto punto no use múltiples núcleos y, por lo tanto, se ejecute ~ 1000x más lento en mi PC.

array: np.random.randint(2, size=shape).astype(dtype)

dtype    shape          %time (average)

float32 (2000, 2000)    62.5 ms
float32 (3000, 3000)    219 ms
float32 (4000, 4000)    328 ms
float32 (10000, 10000)  4.09 s

int8    (2000, 2000)    13 seconds
int8    (3000, 3000)    3min 26s
int8    (4000, 4000)    12min 20s
int8    (10000, 10000)  It didn''t finish in 6 hours

float16 (2000, 2000)    2min 25s
float16 (3000, 3000)    Not tested
float16 (4000, 4000)    Not tested
float16 (10000, 10000)  Not tested

Entiendo que NumPy usa BLAS, que no admite enteros, pero si uso los envoltorios SciPy BLAS, es decir.

import scipy.linalg.blas as blas
a = np.random.randint(2, size=(n, n)).astype(np.int8)
b = blas.sgemm(alpha=1.0, a=a, b=a)

El cálculo es multihilo. Ahora, blas.sgemm ejecuta exactamente con el mismo tiempo que np.dot para float32, pero para los no flotantes convierte todo a float32 y genera flotadores, lo cual es algo que np.dot no hace. (Además, b ahora está en orden F_CONTIGUOUS , que es un problema menor).

Entonces, si quiero hacer la multiplicación de matrices enteras, tengo que hacer uno de los siguientes:

1. Use NumPy, dolorosamente lento, np.dot y esté contento de poder mantener los enteros de 8 bits.
2. Usa el sgemm de SciPy y usa hasta 4x de memoria.
3. Use np.float16 de np.float16 y solo use memoria 2x, con la advertencia de que np.dot es mucho más lento en los arreglos float16 que en los arreglos float32, más que en int8.
4. Encuentre una biblioteca optimizada para la multiplicación de matrices enteras con múltiples subprocesos (en realidad, Mathematica hace esto, pero preferiría una solución de Python), lo ideal sería admitir matrices de 1 bit, aunque las matrices de 8 bits también están bien ... con el objetivo de hacer la multiplicación de matrices sobre el campo finito Z / 2Z, y sé que puedo hacer esto con Sage , que es bastante Pythonic, pero, nuevamente, ¿hay algo estrictamente de Python?)

¿Puedo seguir la opción 4? ¿Existe tal biblioteca?

Descargo de responsabilidad: en realidad estoy ejecutando NumPy + MKL, pero he intentado una prueba similar en NumPy con resultados similares.

Answer 1

" ¿Por qué es más rápido realizar la multiplicación de matrices flotante por flotante en comparación con int por int? " Explica por qué los enteros son tan lentos: en primer lugar, las CPU tienen tuberías de punto flotante de alto rendimiento. En segundo lugar, BLAS no tiene tipo entero.

Solución alternativa: la conversión de matrices a valores float32 consigue grandes aceleraciones. ¿Cómo es la aceleración de 90x en una MacBook Pro 2015? (Usar float64 es la mitad de bueno).

import numpy as np
import time

def timeit(callable):
    start = time.time()
    callable()
    end = time.time()
    return end - start

a = np.random.random_integers(0, 9, size=(1000, 1000)).astype(np.int8)

timeit(lambda: a.dot(a))  # ≈0.9 sec
timeit(lambda: a.astype(np.float32).dot(a.astype(np.float32)).astype(np.int8) )  # ≈0.01 sec

Answer 2

• Opción 5: lanzar una solución personalizada: Particione el producto de matriz en algunos subproductos y realice estos en paralelo. Esto puede ser relativamente fácil de implementar con los módulos estándar de Python. Los subproductos se calculan con numpy.dot , que libera el bloqueo global del intérprete. Por lo tanto, es posible utilizar threads que son relativamente ligeros y pueden acceder a las matrices desde el subproceso principal para la eficiencia de la memoria.

Implementación:

import numpy as np
from numpy.testing import assert_array_equal
import threading
from time import time


def blockshaped(arr, nrows, ncols):
    """
    Return an array of shape (nrows, ncols, n, m) where
    n * nrows, m * ncols = arr.shape.
    This should be a view of the original array.
    """
    h, w = arr.shape
    n, m = h // nrows, w // ncols
    return arr.reshape(nrows, n, ncols, m).swapaxes(1, 2)


def do_dot(a, b, out):
    #np.dot(a, b, out)  # does not work. maybe because out is not C-contiguous?
    out[:] = np.dot(a, b)  # less efficient because the output is stored in a temporary array?


def pardot(a, b, nblocks, mblocks, dot_func=do_dot):
    """
    Return the matrix product a * b.
    The product is split into nblocks * mblocks partitions that are performed
    in parallel threads.
    """
    n_jobs = nblocks * mblocks
    print(''running {} jobs in parallel''.format(n_jobs))

    out = np.empty((a.shape[0], b.shape[1]), dtype=a.dtype)

    out_blocks = blockshaped(out, nblocks, mblocks)
    a_blocks = blockshaped(a, nblocks, 1)
    b_blocks = blockshaped(b, 1, mblocks)

    threads = []
    for i in range(nblocks):
        for j in range(mblocks):
            th = threading.Thread(target=dot_func, 
                                  args=(a_blocks[i, 0, :, :], 
                                        b_blocks[0, j, :, :], 
                                        out_blocks[i, j, :, :]))
            th.start()
            threads.append(th)

    for th in threads:
        th.join()

    return out


if __name__ == ''__main__'':
    a = np.ones((4, 3), dtype=int)
    b = np.arange(18, dtype=int).reshape(3, 6)
    assert_array_equal(pardot(a, b, 2, 2), np.dot(a, b))

    a = np.random.randn(1500, 1500).astype(int)

    start = time()
    pardot(a, a, 2, 4)
    time_par = time() - start
    print(''pardot: {:.2f} seconds taken''.format(time_par))

    start = time()
    np.dot(a, a)
    time_dot = time() - start
    print(''np.dot: {:.2f} seconds taken''.format(time_dot))

Con esta implementación, obtengo una aceleración de aproximadamente x4, que es el número físico de núcleos en mi máquina:

running 8 jobs in parallel
pardot: 5.45 seconds taken
np.dot: 22.30 seconds taken