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

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.

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

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.