Неполная холесская факторизация, очень медленная

Question

Справочная информация: я делаю проект для моего курса по числовой линейной алгебре. Для этого проекта я решил поэкспериментировать с неполной факторизацией Холецкого с арифметикой половинной точности c и использовать результат в качестве предварительного условия для итерационных методов. Сначала я попытался реализовать этот Matlab 2019b (который имеет тип данных с половинной точностью), но он не поддерживает матрицы с половинной точностью разреженные , поэтому мне пришлось использовать полные матрицы. Но арифметика c с половинной точностью в Matlab намного медленнее, и я обнаружил, что для того, чтобы вычислить матрицу 500 x 500, понадобилось около 20 минут (и я хочу подняться до 1000 x 1000). Тем не менее, при одинарной / двойной точности матрица 500 x 500 заняла меньше секунды.

Я подумал, что мне повезло бы с масштабированием на более высокие матрицы, если бы я мог на самом деле воспользоваться разреженностью матрицы. Я вспомнил, что numpy / scipy имеет тип данных float 16, поэтому я решил попробовать реализовать это в python. Поэтому я написал это

from scipy.io import loadmat
def icholesky(a):
    n = a.shape[0]
    for k in tqdm(list(range(n))): 
        a[k,k] = np.sqrt(a[k,k])
        #for i in range(k+1,n):
        #    if (a[i,k] !=0):
        #        a[i,k] = a[i,k]/a[k,k]
        i,_= a[:,k].nonzero()
        if len(i) > 0:
            a[i,k] = a[i,k]/a[k,k]
        for j in range(k+1,n):
            #for i in range(j,n):
            #    if (a[i,j]!=0):
            #        a[i,j] = a[i,j]-a[i,k]*a[j,k]  
            i,_ = a[j:,j].nonzero()
            if len(i) > 0: 
                a[i,j]  = a[i,j] - a[i,k]*a[j,k]     
    return a

bus = loadmat('494_bus.mat') #From University of Florida's Sparse Matrix Collection
A = bus['Problem'][0,0][1]
H = A.copy()
icholesky(H)

Где 'a' будет скудной разреженной матрицей в формате CS C. (Закомментированный код - это просто полностью выписанный алгоритм, не пытающийся воспользоваться преимуществом разреженности). Я обнаружил, что это заняло около 6 минут, что намного быстрее, чем код MATLAB, когда я использую числа с плавающей запятой половинной точности, но все же намного медленнее, чем код Matlab, когда я использую числа с плавающей запятой одинарной / двойной точности (что занимало менее секунды). , хотя MATLAB использует полные матрицы.

Всегда есть вероятность, что я где-то допустил ошибку в своем коде, и я не получаю правильное время выполнения, поэтому я посмотрю это снова. Но мне интересно, если кто-нибудь, более привыкший к scipy / numpy, видит что-то неладное в том, как я решил реализовать приведенный выше код.

У меня есть еще одна теория, почему код python может быть так медленно. Я запускаю это на высокопроизводительном компьютере моей школы, и может случиться так, что matlab настроен на автоматическое использование преимуществ параллелизма, а python - нет. Похоже ли это на разумную гипотезу? Если да, есть ли у вас какие-либо предложения о том, как правильно распараллелить мой алгоритм?

Answer 1

Я нашел способ облегчить проблему. Мне не нужно j, чтобы перебрать все числа от k + 1 до n. Мне просто нужно go через j, где a [j, k] не равно нулю. Что было только что рассчитано, когда я рассчитал i_.

def icholesky(a):
    n = a.shape[0]
    for k in range(n): 
        a[k,k] = np.sqrt(a[k,k])
        i_,_ = a[k+1:,k].nonzero() 
        if len(i_) > 0:
            i_= i_ + (k+1)
            a[i_,k] = a[i_,k]/a[k,k]
        for j in i_: 
            i2_,_ = a[j:n,j].nonzero()
            if len(i2_) > 0:
                i2_ = i2_ + j
                a[i2_,j]  = a[i2_,j] - a[i2_,k]*a[j,k]   


    return a