Почему numpy.linalg.solve () имеет непостоянное время вычисления для массивов одинакового размера? - PullRequest
Новая
       1

Почему numpy.linalg.solve () имеет непостоянное время вычисления для массивов одинакового размера?

0 голосов
/ 18 февраля 2019

Я использовал numpy.linalg.solve(A,B) для решения линейного уравнения.В моем случае: A - это около 10 000x10 000, а B - около 10 000x5.Если я инициализирую A и B случайным образом, используя: 

A = numpy.random.rand(10000,10000)
B = numpy.random.rand(10000,5)

, тогда время вычисления составляет <3 секунды.Тем не менее, в моей программе, которая должна решить это уравнение, вычислительное время постоянно составляет около 14 секунд.Этот код повторяется в цикле, поэтому ускорение почти в 5 раз является большой проблемой.Разве решение для <code>linalg.solve() не должно быть примерно постоянным для массивов постоянного размера?

Обе реализации используют float64.И вполне достаточно оперативной памяти (128 ГБ).Я попытался обновить библиотеки blas на компьютере (Ubuntu 16.04 - numpy установлен с conda), и он показал улучшения в решениях из случайно сгенерированных данных, но не для данных в моей программе.

Если кто-то ищет больше подробностей о коде, это из строки 27 файла .py, найденного в [1].Эта программа выполняет регистрацию облаков точек.

Буду признателен за любые мысли или помощь.

[1] https://github.com/siavashk/pycpd/blob/master/pycpd/deformable_registration.py

Редактировать: Чтобы попытаться сделать это более воспроизводимым, я сгенерировал некоторый код, чтобы попытаться привести меня к подозрительной строке np.linalg.solve (): 

import numpy as np
import time

def gaussian_kernel(Y, beta):
    diff = Y[None,:,:] - Y[:,None,:]
    diff = diff**2
    diff = np.sum(diff, axis=2)
    return np.exp(-diff / (2 * beta**2))

def initialize_sigma2(X, Y):
    diff = X[None,:,:] - Y[:,None,:]
    err  = diff**2
    return np.sum(err) / (X.shape[0] * Y.shape[0] * X.shape[1])

alpha = 0.1
beta = 3

X = np.random.rand(10000,5) * 100
Y = X + X*0.1

N, D = X.shape
M, _ = Y.shape

G = gaussian_kernel(Y, beta)
sigma2 = initialize_sigma2(X,Y)

TY = Y + np.random.rand(10000,5)

P = np.sum((X[None,:,:] - TY[:,None,:])**2, axis=2)
P /= np.sum(P,axis=0)

P1 = np.sum(P, axis=1)
Np = np.sum(P1)

A = np.dot(np.diag(P1), G) + alpha * sigma2 * np.eye(M)
B = np.dot(P, X) - np.dot(np.diag(P1), Y)

%time W = np.linalg.solve(A,B)

Однако этот код не создает временную задержку.Все в этом так же, как в текущем скрипте ... за исключением фактического создания массивов X и Y.Это должны быть облака двух точек, которые приблизительно расположены близко друг к другу в трехмерном пространстве - вот почему я создал одно на основе другого.

...