Comparativo de mercado (python vs. c ++ usando BLAS) e (numpy)

Eu gostaria de escrever um programa que faz uso extensivo das funcionalidades da álgebra linear BLAS e LAPACK. Como o desempenho é um problema, fiz alguns benchmarking e gostaria de saber se a abordagem que usei é legítima.

Tenho, por assim dizer, três concorrentes e quero testar seu desempenho com uma simples multiplicação matriz-matriz. Os competidores são:

1. Numpy, utilizando apenas a funcionalidade do dot.
2. Python, chamando as funcionalidades do BLAS através de um objeto compartilhado.
3. C ++, chamando as funcionalidades do BLAS através de um objeto compartilhado.

Cenário

Implementei uma multiplicação matriz-matriz para diferentes dimensões i. ivai de 5 a 500 com um incremento de 5 e as matrizes m1e m2são configuradas assim:

m1 = numpy.random.rand(i,i).astype(numpy.float32)
m2 = numpy.random.rand(i,i).astype(numpy.float32)

1. Numpy

O código usado é parecido com este:

tNumpy = timeit.Timer("numpy.dot(m1, m2)", "import numpy; from __main__ import m1, m2")
rNumpy.append((i, tNumpy.repeat(20, 1)))

2. Python, chamando BLAS por meio de um objeto compartilhado

Com a função

_blaslib = ctypes.cdll.LoadLibrary("libblas.so")
def Mul(m1, m2, i, r):

    no_trans = c_char("n")
    n = c_int(i)
    one = c_float(1.0)
    zero = c_float(0.0)

    _blaslib.sgemm_(byref(no_trans), byref(no_trans), byref(n), byref(n), byref(n), 
            byref(one), m1.ctypes.data_as(ctypes.c_void_p), byref(n), 
            m2.ctypes.data_as(ctypes.c_void_p), byref(n), byref(zero), 
            r.ctypes.data_as(ctypes.c_void_p), byref(n))

o código de teste é semelhante a este:

r = numpy.zeros((i,i), numpy.float32)
tBlas = timeit.Timer("Mul(m1, m2, i, r)", "import numpy; from __main__ import i, m1, m2, r, Mul")
rBlas.append((i, tBlas.repeat(20, 1)))

3. c ++, chamando BLAS por meio de um objeto compartilhado

Agora o código c ++ naturalmente é um pouco mais longo, então reduzo as informações ao mínimo.
Eu carrego a função com

void* handle = dlopen("libblas.so", RTLD_LAZY);
void* Func = dlsym(handle, "sgemm_");

Eu medi o tempo gettimeofdayassim:

gettimeofday(&start, NULL);
f(&no_trans, &no_trans, &dim, &dim, &dim, &one, A, &dim, B, &dim, &zero, Return, &dim);
gettimeofday(&end, NULL);
dTimes[j] = CalcTime(start, end);

onde jestá um loop rodando 20 vezes. Eu calculo o tempo que passou com

double CalcTime(timeval start, timeval end)
{
double factor = 1000000;
return (((double)end.tv_sec) * factor + ((double)end.tv_usec) - (((double)start.tv_sec) * factor + ((double)start.tv_usec))) / factor;
}

Resultados

O resultado é mostrado no gráfico abaixo:

Questões

1. Você acha que minha abordagem é justa ou há alguma sobrecarga desnecessária que posso evitar?
2. Você esperaria que o resultado mostrasse uma discrepância tão grande entre a abordagem c ++ e python? Ambos estão usando objetos compartilhados para seus cálculos.
3. Como prefiro usar python para meu programa, o que posso fazer para aumentar o desempenho ao chamar rotinas BLAS ou LAPACK?

Baixar

O benchmark completo pode ser baixado aqui . (JF Sebastian tornou esse link possível ^^)

Eu executei seu benchmark . Não há diferença entre C ++ e numpy na minha máquina:

Você acha que minha abordagem é justa ou há alguma sobrecarga desnecessária que posso evitar?

Parece justo porque não há diferença nos resultados.

Você esperaria que o resultado mostrasse uma discrepância tão grande entre a abordagem c ++ e python? Ambos estão usando objetos compartilhados para seus cálculos.

Não.

Como prefiro usar python para meu programa, o que posso fazer para aumentar o desempenho ao chamar rotinas BLAS ou LAPACK?

Certifique-se de que numpy usa uma versão otimizada das bibliotecas BLAS / LAPACK em seu sistema.

ATUALIZAÇÃO (30.07.2014):

Eu executo novamente o benchmark em nosso novo HPC. Tanto o hardware quanto a pilha de software mudaram da configuração na resposta original.

Coloquei os resultados em uma planilha do Google (contém também os resultados da resposta original).

Hardware

Nosso HPC tem dois nós diferentes, um com CPUs Intel Sandy Bridge e outro com as CPUs Ivy Bridge mais recentes:

Sandy (MKL, OpenBLAS, ATLAS):

• CPU : 2 x 16 Intel (R) Xeon (R) E2560 Sandy Bridge a 2,00 GHz (16 núcleos)
• RAM : 64 GB

Ivy (MKL, OpenBLAS, ATLAS):

• CPU : 2 x 20 Intel (R) Xeon (R) E2680 V2 Ivy Bridge @ 2,80 GHz (20 núcleos, com HT = 40 núcleos)
• RAM : 256 GB

Programas

A pilha de software é para ambos os nós o sam. Em vez de GotoBLAS2 , OpenBLAS é usado e também há um ATLAS BLAS multithread que é definido para 8 threads (codificado).

• SO : Suse
• Compilador Intel : ictce-5.3.0
• Numpy: 1.8.0
• OpenBLAS: 0.2.6
• ATLAS :: 3.8.4

Referência de produto pontual

O código de referência é o mesmo que abaixo. No entanto, para as novas máquinas, também executei o benchmark para os tamanhos de matriz 5000 e 8000 .
A tabela abaixo inclui os resultados do benchmark da resposta original (renomeado: MKL -> Nehalem MKL, Netlib Blas -> Nehalem Netlib BLAS, etc)

Desempenho de thread único:

Desempenho multi thread (8 threads):

Threads vs tamanho da matriz (Ivy Bridge MKL) :

Suite Benchmark

Desempenho de thread único:

Desempenho multiencadeado (8 fios):

Conclusão

Os novos resultados do benchmark são semelhantes aos da resposta original. OpenBLAS e MKL funcionam no mesmo nível, com exceção do teste de autovalor . Os Eigenvalue executa o teste única razoavelmente bem no OpenBLAS em modo de rosca única . No modo multithread, o desempenho é pior.

O "gráfico de tamanho da matriz vs threads" também mostra que, embora o MKL, assim como o OpenBLAS, geralmente escalem bem com o número de núcleos / threads, isso depende do tamanho da matriz. Para matrizes pequenas, adicionar mais núcleos não melhora muito o desempenho.

Há também um aumento de desempenho de aproximadamente 30% de Sandy Bridge para Ivy Bridge, o que pode ser devido a uma freqüência maior (+ 0,8 Ghz) e / ou melhor arquitetura.

Resposta Original (04.10.2011):

Algum tempo atrás eu tive que otimizar alguns cálculos / algoritmos de álgebra linear que foram escritos em python usando numpy e BLAS, então eu comparei / testei diferentes configurações numpy / BLAS.

Eu testei especificamente:

• Numpy com ATLAS
• Numpy com GotoBlas2 (1.13)
• Numpy com MKL (11.1 / 073)
• Numpy com Accelerate Framework (Mac OS X)

Eu executei dois benchmarks diferentes:

1. produto escalar simples de matrizes com tamanhos diferentes
2. Suite de referência que pode ser encontrada aqui .

Aqui estão meus resultados:

Maquinas

Linux (MKL, ATLAS, No-MKL, GotoBlas2):

• SO : Ubuntu Lucid 10.4 64 bits.
• CPU : 2 x 4 Intel (R) Xeon (R) E5504 @ 2,00 GHz (8 núcleos)
• RAM : 24 GB
• Compilador Intel : 11.1 / 073
• Scipy : 0,8
• Numpy : 1,5

Mac Book Pro (Accelerate Framework):

• SO : Mac OS X Snow Leopard (10.6)
• CPU : 1 Intel Core 2 Duo 2,93 Ghz (2 núcleos)
• RAM : 4 GB
• Scipy : 0,7
• Numpy : 1,3

Servidor Mac (Accelerate Framework):

• SO : Servidor Mac OS X Snow Leopard (10.6)
• CPU : 4 X Intel (R) Xeon (R) E5520 a 2,26 Ghz (8 núcleos)
• RAM : 4 GB
• Scipy : 0,8
• Numpy : 1.5.1

Referência de produto pontual

Código :

import numpy as np
a = np.random.random_sample((size,size))
b = np.random.random_sample((size,size))
%timeit np.dot(a,b)

Resultados :

    Sistema | size = 1000 | tamanho = 2000 | size = 3000 |
netlib BLAS | 1350 ms | 10900 ms | 39200 ms |    
ATLAS (1 CPU) | 314 ms | 2560 ms | 8700 ms |     
MKL (1 CPU) | 268 ms | 2110 ms | 7120 ms |
MKL (2 CPUs) | - | - | 3660 ms |
MKL (8 CPUs) | 39 ms | 319 ms | 1000 ms |
GotoBlas2 (1 CPU) | 266 ms | 2100 ms | 7280 ms |
GotoBlas2 (2 CPUs) | 139 ms | 1009 ms | 3690 ms |
GotoBlas2 (8 CPUs) | 54 ms | 389 ms | 1250 ms |
Mac OS X (1 CPU) | 143 ms | 1060 ms | 3605 ms |
Servidor Mac (1 CPU) | 92 ms | 714 ms | 2130 ms |

Suite Benchmark

Código :
para obter informações adicionais sobre o pacote de benchmarks, clique aqui .

Resultados :

    Sistema | valores próprios | svd | det | inv | ponto |
netlib BLAS | 1688 ms | 13102 ms | 438 ms | 2155 ms | 3522 ms |
ATLAS (1 CPU) | 1210 ms | 5897 ms | 170 ms | 560 ms | 893 ms |
MKL (1 CPU) | 691 ms | 4475 ms | 141 ms | 450 ms | 736 ms |
MKL (2 CPUs) | 552 ms | 2718 ms | 96 ms | 267 ms | 423 ms |
MKL (8 CPUs) | 525 ms | 1679 ms | 60 ms | 137 ms | 197 ms |  
GotoBlas2 (1 CPU) | 2124 ms | 4636 ms | 147 ms | 456 ms | 743 ms |
GotoBlas2 (2 CPUs) | 1560 ms | 3278 ms | 116 ms | 295 ms | 460 ms |
GotoBlas2 (8 CPUs) | 741 ms | 2914 ms | 82 ms | 262 ms | 192 ms |
Mac OS X (1 CPU) | 948 ms | 4339 ms | 151 ms | 318 ms | 566 ms |
Servidor Mac (1 CPU) | 1033 ms | 3645 ms | 99 ms | 232 ms | 342 ms |

Instalação

A instalação do MKL incluiu a instalação do Intel Compiler Suite completo, que é bastante simples. No entanto, por causa de alguns bugs / problemas, configurar e compilar numpy com suporte MKL foi um pouco incômodo.

GotoBlas2 é um pequeno pacote que pode ser facilmente compilado como uma biblioteca compartilhada. No entanto, por causa de um bug, você deve recriar a biblioteca compartilhada após compilá-la para usá-la com o numpy.
Além dessa construção para plataforma de múltiplos alvos, não funcionou por algum motivo. Então, eu tive que criar um arquivo .so para cada plataforma para a qual quero ter um arquivo libgoto2.so otimizado .

Se você instalar o numpy do repositório do Ubuntu, ele irá instalar e configurar automaticamente o numpy para usar o ATLAS . A instalação do ATLAS a partir da fonte pode demorar algum tempo e requer algumas etapas adicionais (fortran, etc).

Se você instalar o numpy em uma máquina Mac OS X com portas Fink ou Mac, ele configurará o numpy para usar ATLAS ou o Accelerate Framework da Apple . Você pode verificar executando ldd no arquivo numpy.core._dotblas ou chamando numpy.show_config () .

Conclusões

MKL tem melhor desempenho, seguido de perto por GotoBlas2 .
No teste de autovalor , GotoBlas2 tem um desempenho surpreendentemente pior do que o esperado. Não sei por que esse é o caso.
O Accelerate Framework da Apple tem um desempenho muito bom, especialmente no modo single threaded (em comparação com as outras implementações BLAS).

Ambos GotoBlas2 e MKL escalam muito bem com o número de threads. Portanto, se você tiver que lidar com grandes matrizes, executá-lo em vários threads ajudará muito.

Em qualquer caso, não use a implementação padrão do netlib blas porque é muito lenta para qualquer trabalho computacional sério.

Em nosso cluster, também instalei o ACML da AMD e o desempenho foi semelhante ao MKL e GotoBlas2 . Eu não tenho nenhum número difícil.

Eu pessoalmente recomendaria usar GotoBlas2 porque é mais fácil de instalar e é grátis.

Se você quiser codificar em C ++ / C também verifique o Eigen3, que supostamente supera o MKL / GotoBlas2 em alguns casos e também é muito fácil de usar.

Aqui está outra referência (no Linux, basta digitar make): http://dl.dropbox.com/u/5453551/blas_call_benchmark.zip

http://dl.dropbox.com/u/5453551/blas_call_benchmark.png

Não vejo essencialmente nenhuma diferença entre os diferentes métodos para grandes matrizes, entre Numpy, Ctypes e Fortran. (Fortran em vez de C ++ --- e se isso for importante, seu benchmark provavelmente está quebrado.)

Sua CalcTimefunção em C ++ parece ter um erro de sinal. ... + ((double)start.tv_usec))deve ser em vez disso ... - ((double)start.tv_usec)). Talvez seu benchmark também tenha outros bugs, por exemplo, comparar entre diferentes bibliotecas BLAS ou diferentes configurações BLAS, como número de threads, ou entre tempo real e tempo de CPU?

EDIT : falha ao contar as chaves na CalcTimefunção - está OK.

Como uma diretriz: se você fizer um benchmark, sempre poste todo o código em algum lugar. Comentar sobre benchmarks, especialmente quando surpreendente, sem ter o código completo geralmente não é produtivo.

Para descobrir contra qual BLAS Numpy está vinculado, faça:

$ python
Python 2.7.2+ (padrão, 16 de agosto de 2011, 07:24:41) 
[GCC 4.6.1] no linux2
Digite "ajuda", "copyright", "créditos" ou "licença" para obter mais informações.
>>> import numpy.core._dotblas
>>> numpy.core._dotblas .__ file__
'/usr/lib/pymodules/python2.7/numpy/core/_dotblas.so'
>>> 
$ ldd /usr/lib/pymodules/python2.7/numpy/core/_dotblas.so
    linux-vdso.so.1 => (0x00007fff5ebff000)
    libblas.so.3gf => /usr/lib/libblas.so.3gf (0x00007fbe618b3000)
    libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007fbe61514000)

ATUALIZAÇÃO : Se você não pode importar numpy.core._dotblas, seu Numpy está usando sua cópia interna do BLAS, que é mais lenta, e não deve ser usada em computação de desempenho! A resposta de @Woltan abaixo indica que esta é a explicação para a diferença que ele vê em Numpy vs. Ctypes + BLAS.

Para corrigir a situação, você precisa de ATLAS ou MKL --- verifique estas instruções: http://scipy.org/Installing_SciPy/Linux A maioria das distribuições Linux vem com ATLAS, então a melhor opção é instalar o libatlas-devpacote (o nome pode variar) .

Dado o rigor que você mostrou com sua análise, estou surpreso com os resultados até agora. Eu coloquei isso como uma 'resposta', mas apenas porque é muito longo para um comentário e fornece uma possibilidade (embora eu espere que você tenha considerado).

Eu teria pensado que a abordagem numpy / python não acrescentaria muito overhead para uma matriz de complexidade razoável, já que à medida que a complexidade aumenta, a proporção em que o python participa deve ser pequena. Estou mais interessado nos resultados do lado direito do gráfico, mas a discrepância de ordens de magnitude mostrada lá seria preocupante.

Eu me pergunto se você está usando os melhores algoritmos que o entorpecimento pode aproveitar. Do guia de compilação para Linux:

"Build FFTW (3.1.2): Versões SciPy> = 0.7 e Numpy> = 1.2: Por causa dos problemas de licença, configuração e manutenção, o suporte para FFTW foi removido nas versões de SciPy> = 0.7 e NumPy> = 1.2. Agora usa uma versão integrada do fftpack. Existem algumas maneiras de aproveitar a velocidade do FFTW, se necessário, para sua análise. Faça downgrade para uma versão Numpy / Scipy que inclua suporte. Instale ou crie seu próprio wrapper do FFTW. Consulte http: //developer.berlios.de/projects/pyfftw/ como um exemplo não aprovado. "

Você compilou numpy com mkl? ( http://software.intel.com/en-us/articles/intel-mkl/ ). Se você estiver executando no Linux, as instruções para compilar numpy com mkl estão aqui: http://www.scipy.org/Installing_SciPy/Linux#head-7ce43956a69ec51c6f2cedd894a4715d5bfff974 (apesar de url). A parte principal é:

[mkl]
library_dirs = /opt/intel/composer_xe_2011_sp1.6.233/mkl/lib/intel64
include_dirs = /opt/intel/composer_xe_2011_sp1.6.233/mkl/include
mkl_libs = mkl_intel_lp64,mkl_intel_thread,mkl_core 

Se você estiver no Windows, pode obter um binário compilado com mkl (e também obter o pyfftw e muitos outros algoritmos relacionados) em: http://www.lfd.uci.edu/~gohlke/pythonlibs/ , com um dívida de gratidão para com Christoph Gohlke no Laboratory for Fluorescence Dynamics, UC Irvine.

Advertência, em qualquer caso, há muitos problemas de licenciamento e assim por diante para estar ciente, mas a página Intel explica esses. Novamente, imagino que você tenha considerado isso, mas se você atender aos requisitos de licenciamento (o que no Linux é muito fácil de fazer), isso aceleraria muito a parte entorpecente em relação ao uso de uma construção automática simples, sem nem mesmo FFTW. Terei interesse em seguir este tópico e ver o que os outros pensam. Apesar de tudo, excelente rigor e excelente questão. Obrigado por postar.