Por que o np.dot é impreciso? (matrizes n-dim)

15

Suponha que tomemos np.dotduas 'float32'matrizes 2D:

res = np.dot(a, b)   # see CASE 1
print(list(res[0]))  # list shows more digits
[-0.90448684, -1.1708503, 0.907136, 3.5594249, 1.1374011, -1.3826287]

Números. Exceto, eles podem mudar:


CASO 1 : fatiaa

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0])) # full shape: (i, 6)
[-0.9044868,  -1.1708502, 0.90713596, 3.5594249, 1.1374012, -1.3826287]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]

Os resultados diferem, mesmo que a fatia impressa derive exatamente dos mesmos números multiplicados.


CASO 2 : achatar a, dar uma versão 1D de b, em seguida, fatia a:

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(1, 6).astype('float32')

for i in range(1, len(a)):
    a_flat = np.expand_dims(a[:i].flatten(), -1) # keep 2D
    print(list(np.dot(a_flat, b)[0])) # full shape: (i*6, 6)
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]

CASO 3 : controle mais forte; defina todos os itens não envolvidos como zero : adicione a[1:] = 0ao código do CASO 1. Resultado: discrepâncias persistem.


CASO 4 : verifique índices diferentes de [0]; como para [0], os resultados começam a estabilizar um número fixo de ampliações de matriz a partir do ponto de criação. Resultado

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for j in range(len(a) - 2):
    for i in range(1, len(a)):
        res = np.dot(a[:i], b)
        try:    print(list(res[j]))
        except: pass
    print()

Portanto, para o caso 2D * 2D, os resultados diferem - mas são consistentes para 1D * 1D. De algumas das minhas leituras, isso parece resultar de 1D-1D usando adição simples, enquanto 2D-2D usa adição 'mais sofisticada', que aumenta o desempenho, que pode ser menos precisa (por exemplo, a adição pareada faz o oposto). No entanto, não consigo entender por que as discrepâncias desaparecem no caso em que 1 aé fatiado além de um 'limiar' definido; o maior aeb mais tarde esse limite parece estar, mas sempre existe.

Tudo dito: por que é np.dotimpreciso (e inconsistente) para matrizes ND-ND? Git Relevante


Informações adicionais :

  • Meio Ambiente : operacional Win-10, Python 3.7.4, Spyder 3.3.6 IDE, Anaconda 3.0 2019/10
  • CPU : i7-7700HQ 2,8 GHz
  • Numpy v1.16.5

Possível biblioteca culpada : Numpy MKL - também bibliotecas BLASS; graças a Bi Rico por notar


Código do teste de estresse : como observado, as discrepâncias exacerbam na frequência com matrizes maiores; se acima não for reproduzível, abaixo deve ser (se não, tente dims maiores). Minha saída

np.random.seed(1)
a = (0.01*np.random.randn(9, 9999)).astype('float32') # first multiply then type-cast
b = (0.01*np.random.randn(9999, 6)).astype('float32') # *0.01 to bound mults to < 1

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0]))

Gravidade do problema : as discrepâncias mostradas são 'pequenas', mas não mais quando operam em uma rede neural com bilhões de números multiplicados por alguns segundos e trilhões por todo o tempo de execução; a precisão do modelo relatado difere por 10 inteiros de porcentagens, de acordo com este segmento .

Abaixo está um gif de matrizes resultantes da alimentação para um modelo, basicamente a[0], w / len(a)==1vs len(a)==32.:


OUTRAS PLATAFORMAS resultados, de acordo e com agradecimentos a Paul teste de :

Caso 1 reproduzido (parcialmente) :

  • VM do Google Colab - Intel Xeon 2.3 G-Hz - Jupyter - Python 3.6.8
  • Win-10 Pro Docker Desktop - Intel i7-8700K - jupyter / scipy-notebook - Python 3.7.3
  • Ubuntu 18.04.2 LTS + Docker - AMD FX-8150 - jupyter / scipy-notebook - Python 3.7.3

Nota : estes geram um erro muito menor do que o mostrado acima; duas entradas na primeira linha são desativadas por 1 no dígito menos significativo das entradas correspondentes nas outras linhas.

Caso 1 não reproduzido :

  • Ubuntu 18.04.3 LTS - Intel i7-8700K - IPython 5.5.0 - Python 2.7.15+ e 3.6.8 (2 testes)
  • Ubuntu 18.04.3 LTS - Intel i5-3320M - IPython 5.5.0 - Python 2.7.15+
  • Ubuntu 18.04.2 LTS - AMD FX-8150 - IPython 5.5.0 - Python 2.7.15rc1

Notas :

  • O ambiente de notebook e jupiter da Colab vinculado mostra uma discrepância muito menor (e apenas nas duas primeiras linhas) do que é observado no meu sistema. Além disso, o Caso 2 nunca (ainda) mostrou imprecisão.
  • Dentro desta amostra muito limitada, o atual ambiente Jupyter (Dockerized) é mais suscetível que o ambiente IPython.
  • np.show_config()muito tempo para postar, mas em resumo: os envs do IPython são baseados em BLAS / LAPACK; O Colab é baseado no OpenBLAS. Nos envs do IPython Linux, as bibliotecas BLAS são instaladas pelo sistema - no Jupyter e no Colab, elas vêm de / opt / conda / lib

ATUALIZAÇÃO : a resposta aceita é precisa, mas ampla e incompleta. A questão permanece em aberto para qualquer um que possa explicar o comportamento no nível do código - ou seja, um algoritmo exato usado por np.dot, e como ele explica 'inconsistências consistentes' observadas nos resultados acima (veja também comentários). Aqui estão algumas implementações diretas além da minha decifração: sdot.c - arraytypes.c.src

OverLordGoldDragon
fonte
Comentários não são para discussão prolongada; esta conversa foi movida para o bate-papo .
Samuel Liew
Os algoritmos gerais para ndarraysgeralmente desconsideram a perda de precisão numérica. Porque a simplicidade que reduce-sumao longo de cada eixo, a ordem das operações pode não ser o ideal um ... Note que se você mente erro de precisão assim como você pode usarfloat64
Vitor SRG
Talvez eu não tenha tempo para revisar amanhã, então conceda a recompensa agora.
Paul
@Paul Seria concedido automaticamente à resposta mais votada de qualquer maneira - mas tudo bem, obrigado por notificar
OverLordGoldDragon

Respostas:

7

Parece imprecisão numérica inevitável. Como explicado aqui , o NumPy usa um método BLAS altamente otimizado e cuidadosamente ajustado para multiplicação de matrizes . Isso significa que provavelmente a sequência de operações (soma e produtos) seguida para multiplicar 2 matrizes muda quando o tamanho da matriz muda.

Tentando ser mais claro, sabemos que, matematicamente , cada elemento da matriz resultante pode ser calculado como o produto escalar de dois vetores (sequências de números de comprimento igual). Mas não é assim que o NumPy calcula um elemento da matriz resultante. De fato, existem algoritmos mais eficientes, mas complexos, como o algoritmo Strassen , que obtêm o mesmo resultado sem calcular diretamente o produto de ponto da coluna de linha.

Ao usar esses algoritmos, mesmo que o elemento C ij de uma matriz resultante C = AB seja matematicamente definido como o produto escalar da i-ésima fileira de um com o j-th coluna de B , se multiplicar uma matriz A2 ter o mesma i-ésima linha que A com uma matriz B2 com a mesma j-ésima coluna de B , o elemento C2 ij será realmente calculado seguindo uma sequência diferente de operações (que depende de todo A2 e B2 matrizes), possivelmente levando a diferentes erros numéricos.

É por isso que, mesmo que matematicamente C ij = C2 ij (como em seu caso 1), a diferente sequência de operações seguido pelo algoritmo nos cálculos (devido à mudança no tamanho da matriz) leva a diferentes erros numéricos. O erro numérico explica também os resultados ligeiramente diferentes, dependendo do ambiente e o fato de que, em alguns casos, em alguns ambientes, o erro numérico pode estar ausente.

mmj
fonte
2
Obrigado pelo link, ele parece conter informações relevantes - sua resposta, no entanto, pode ser mais detalhada, pois até agora é uma paráfrase dos comentários abaixo da pergunta. Por exemplo, o SO vinculado mostra Ccódigo direto e fornece explicações em nível de algoritmo, por isso está indo na direção certa.
OverLordGoldDragon
Também não é "inevitável", como mostrado na parte inferior da pergunta - e a extensão da imprecisão varia entre os ambientes, o que permanece inexplicável
OverLordGoldDragon
11
@OverLordGoldDragon: (1) Um exemplo trivial com a adição: pegue o número n, pegue o número de kforma que fique abaixo da precisão do kúltimo dígito da mantissa. Para carros alegóricos nativos do Python, n = 1.0e k = 1e-16funciona. Agora vamos ks = [k] * 100. Veja que sum([n] + ks) == n, enquanto sum(ks + [n]) > nisso é importante, a ordem de soma. (2) As CPUs modernas têm várias unidades para executar operações de ponto flutuante (FP) em paralelo, e a ordem na qual a + b + c + dé computada em uma CPU não é definida, mesmo se o comando a + bvier antes c + dno código da máquina.
9000
11
@OverLordGoldDragon Você deve estar ciente de que a maioria dos números com os quais você solicita que seu programa lide não pode ser representada exatamente por um ponto flutuante. Tente format(0.01, '.30f'). Se mesmo um número simples como 0.01não puder ser representado exatamente por um ponto flutuante NumPy, não há necessidade de conhecer os detalhes profundos do algoritmo de multiplicação de matrizes NumPy para entender o ponto da minha resposta; isto é, matrizes iniciais diferentes levam a diferentes seqüências de operações , para que resultados matematicamente iguais possam diferir em uma pequena quantidade devido a erros numéricos.
MMJ
2
@OverLordGoldDragon re: magia negra. Há um artigo que é leitura obrigatória em alguns MOOCs da CS. Minha recordação não é tão grande, mas eu acho que é isso: itu.dk/~sestoft/bachelor/IEEE754_article.pdf
Paul