Полная диагонализация NumPy против SciPy

Я должен использовать процедуры диагонализации, чтобы получить все собственные пары эрмитовой комплексной матрицы. Я немного ограничен производительностью, так как мне нужно повторять операцию тысячи раз, а мои матрицы примерно 8000x8000. Я создал небольшое сравнение между подпрограммами NumPy и SciPy для диагонализации эрмитовых матриц и получил эти времена на машине с 6 физическими ядрами:

enter image description here

Я наблюдаю, что для матриц 8000x8000 это масштабируется до ~ 0,8 минуты, и мне нужно повторить процесс 50000 раз. Есть ли что-то, что мне здесь не хватает, или это на самом деле время выступления? В целом, все это выглядит очень медленно, особенно если это нужно повторить несколько раз. Фактически, на 30-ядерной машине я наблюдаю небольшой прирост производительности. Я использую Python3.8 в дистрибутиве Anaconda, поэтому он связан с MKL.

Вот пример кода

import numpy as np
import scipy.linalg
import matplotlib.pyplot as pyt
from time import time


t_ls = []
d_ls = np.array([100, 500, 1000,2000,4000])

for N in d_ls:

    A =np.random.rand( N,N ) + 1j*np.random.rand( N,N )    
    A = 0.5*( A + np.conj(A.T) ) 
       
    ts = time()    
    evals, evecs = np.linalg.eigh( A )
    t_np  = time()-ts
   

    ts = time()    
    evals2, evecs2 = scipy.linalg.eigh( A )
    t_sp  = time()-ts
    
    t_ls.append(np.array([t_np, t_sp]))
    
t_ls = np.array(t_ls)

pyt.plot( d_ls, t_ls[:,0], marker='s' )
pyt.plot( d_ls, t_ls[:,1], marker='^')
pyt.xlabel("N")
pyt.ylabel("time(secs)")
pyt.legend(["NumPy", "SciPy"])
pyt.show()

ИСПОЛЬЗОВАНИЕ СВД И МП ПАРАЛЕЛЛИЗАЦИЯ

Просматривая некоторые комментарии в посте, я попробовал SVD матрицы и многопроцессорности. Во всех случаях я все еще вижу, что сериализованный подход с NumPy eigh является наиболее эффективным; вот код:

import numpy as np
import scipy.linalg
import matplotlib.pyplot as pyt
from time import time

import psutil


def f_mp_pool(*args):
     
    N = args[0]
        
    A =np.random.rand( N,N ) + 1j*np.random.rand( N,N )    
    A = 0.5*( A + np.conj(A.T) ) 
    
    evals, evecs = np.linalg.eigh(A)
    
    return evals, evecs


nreps = 100
N     = 700


ts = time()

for n in range(nreps):
    
    A =np.random.rand( N,N ) + 1j*np.random.rand( N,N )    
    A = 0.5*( A + np.conj(A.T) ) 
    
    res = np.linalg.eigh(A)
    

print("serialized:", time()-ts)

#use svd

import scipy.linalg

ts = time()
for n in range(nreps):
    res = scipy.linalg.svd( A, full_matrices=True, check_finite=False  ) 
    
print("SVD:", time()-ts)    

import multiprocessing as mp

nproc   = psutil.cpu_count(logical=False)-1
mp_pool = mp.Pool(processes=nproc)

args_ls = [ (N,) for n in range(nreps) ]

ts = time()
res = mp_pool.starmap( f_mp_pool, args_ls )
print("parallel:", time()-ts)

Неудивительно, что вы видите небольшое ускорение независимо от количества ядер, поскольку вы выполняете все процессы последовательно, то есть один за другим. Чтобы действительно использовать несколько ядер, взгляните на многопроцессорность.

DrBwts 17.05.2022 10:25

Но я всегда думал, что многопроцессорность чего-то, что использует все потоки (например, из-за MKL), не принесет никаких преимуществ. Итак, если у меня 30 ядер, и я распараллеливаю, как вы предлагаете, устанавливая количество потоков равным 2 для каждого запущенного процесса, не будет ли это еще медленнее?

Zarathustra 17.05.2022 10:31

Это похоже на хороший пример использования Numba, Cython или даже vmap JAX.

keepAlive 17.05.2022 10:38

Как Numba помогает с диагонализацией? Разве это не просто обертка всей процедуры NumPy? См. stackoverflow.com/questions/57603502/…

Zarathustra 17.05.2022 10:46

Почему бы просто не запустить SVD вместо ED?

onyambu 17.05.2022 11:39

Я пробовал использовать SVD, но это не работает быстрее; дополнительно нужно взять квадрат для собственных значений и т.д...

Zarathustra 17.05.2022 11:53

хм, я продолжаю получать ошибку времени выполнения, когда ваш код входит в часть multiprocessing

DrBwts 17.05.2022 18:32

Можно уточнить какая ошибка? Минимально воспроизводимый код отлично работает на моей машине.

Zarathustra 17.05.2022 21:05
Формы c голосовым вводом в React с помощью Speechly
Формы c голосовым вводом в React с помощью Speechly
Пытались ли вы когда-нибудь заполнить веб-форму в области электронной коммерции, которая требует много кликов и выбора? Вас попросят заполнить дату,...
Стилизация и валидация html-формы без использования JavaScript (только HTML/CSS)
Стилизация и валидация html-формы без использования JavaScript (только HTML/CSS)
Будучи разработчиком веб-приложений, легко впасть в заблуждение, считая, что приложение без JavaScript не имеет права на жизнь. Нам становится удобно...
Flatpickr: простой модуль календаря для вашего приложения на React
Flatpickr: простой модуль календаря для вашего приложения на React
Если вы ищете пакет для быстрой интеграции календаря с выбором даты в ваше приложения, то библиотека Flatpickr отлично справится с этой задачей....
В чем разница между Promise и Observable?
В чем разница между Promise и Observable?
Разберитесь в этом вопросе, и вы значительно повысите уровень своей компетенции.
Что такое cURL в PHP? Встроенные функции и пример GET запроса
Что такое cURL в PHP? Встроенные функции и пример GET запроса
Клиент для URL-адресов, cURL, позволяет взаимодействовать с множеством различных серверов по множеству различных протоколов с синтаксисом URL.
Четыре эффективных способа центрирования блочных элементов в CSS
Четыре эффективных способа центрирования блочных элементов в CSS
У каждого из нас бывали случаи, когда нам нужно отцентрировать блочный элемент, но мы не знаем, как это сделать. Даже если мы реализуем какой-то...
0
8
34
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Питорч будет работать быстрее, и если у вас есть GPU, он также может воспользоваться этим преимуществом, но не так сильно, потому что итерация QR не подходит для параллельных вычислений. У меня есть потенциальное решение для ускорения этой части на графических процессорах, но я так и не реализовал его.

import numpy as np
import scipy.linalg
import torch
import matplotlib.pyplot as plt
from time import time

t_ls = []
d_ls = np.array([100, 500, 1000,2000,4000])

for N in d_ls:

    A =np.random.rand( N,N ) + 1j*np.random.rand( N,N )    
    A = 0.5*( A + np.conj(A.T) ) 

#  skipping numpy, it is slow here, you may put it back if you want
#     ts = time()    
#     evals, evecs = np.linalg.eigh( A )
#     t_np  = time()-ts
   

    ts = time()    
    evals2, evecs2 = scipy.linalg.eigh( A )
    t_sp  = time()-ts
    
    # When using CPU torch will use intra operation
    # parallelism, so if you care about latency
    # this is better than using multiprocessing
    A_cpu = torch.as_tensor(A)
    ts = time()
    evals3, evecs3 = torch.linalg.eigh(A_cpu)
    t_cpu = time() - ts;
    if torch.cuda.is_available():
        # Using GPU will give a significant speedup for some
        # operations, I guess the 
        A_gpu = A_cpu.to('cuda')
        torch.cuda.synchronize()
        ts = time()
        evals4, evecs4 = torch.linalg.eigh(A_gpu)
        torch.cuda.synchronize()
        t_gpu = time() - ts;
    else:
        t_gpu = np.nan #if you don't have GPU let's skip this part
    t_ls.append(np.array([np.nan, t_sp, t_cpu, t_gpu]))
    print(t_ls[-1])

t_ls = np.array(t_ls)

plt.plot( d_ls, t_ls[:,0], marker='s' )
plt.plot( d_ls, t_ls[:,1], marker='^')
plt.plot( d_ls, t_ls[:,2], marker='+')
plt.plot( d_ls, t_ls[:,3], marker='d')
plt.xlabel("N")
plt.ylabel("time(secs)")
plt.legend(["NumPy", "SciPy", "PyTorch CPU", "PyTorch GPU"])

Мой сюжет enter image description here

Это выглядит как большое улучшение. Я недавно читал о Torch и действительно только что попробовал обычный подход, но GPU значительно повышает производительность! Кое-что, что я также отметил в многопроцессорном распараллеливании, заключается в том, что установка числа потоков MKL равным 1, похоже, делает работу намного быстрее; каждый процесс затем ограничивается одним процессором.

Zarathustra 17.05.2022 21:05

Если вы хотите одновременно вычислять несколько матриц одинакового размера, вы можете рассмотреть возможность их укладки, pytorch работает с пакетом матриц. Я знаю, что для низких размеров это не очень помогает, но, может быть, для вашего варианта использования.

Bob 17.05.2022 22:47

Другие вопросы по теме