← Materiały bloku: Kwantowe wyżarzanie kombinatorycznych problemów optymalizacyjnych  ·  podgląd notatnika 06_GPU.ipynb  ·  pobierz plik .ipynb

Wykorzystanie procesorów graficznych (GPU) w algorytmach heurystycznych

Wstęp - jak używać GPU

Podstawową zaletą używania GPU jest duża możliwość równoległego wykonywania obliczeń, co pozwala znacząco przyspieszyć operacje wymagające intensywnej mocy obliczeniowej. W praktyce istnieją dwa główne sposoby wykorzystania GPU: można sięgnąć po gotowe implementacje dostępne w popularnych bibliotekach i frameworkach (np. TensorFlow, PyTorch czy CuPy), które automatycznie obsługują akcelerację na karcie graficznej, lub napisać własne jądro CUDA, co daje pełną kontrolę nad strukturą obliczeń i umożliwia optymalizację pod kątem konkretnego zadania.

Pisanie własnych jąder CUDA wymaga pewnej wprawy, ale pozwala na uzyskanie maksymalnej wydajności w niestandardowych zastosowaniach. Pardzo dobry poradnik jak to robić mozna znaleść na stronie CUDA C++ Programming Guide.

Symulowana bifurkacja

# Balistyczna symulowana bifurkacja na gpu naiwna - zamieniamy np na cp

from setup import setup_paths
setup_paths()


import cupy as cp
from math import sqrt
from typing import Optional
from tqdm import tqdm

from funkcje_pomocnicze_gpu import calculate_energy_gpu


def wall_gpu(x: cp.ndarray, y: cp.ndarray):
    mask = cp.abs(x) > 1
    x[mask] = cp.sign(x[mask])
    y[mask] = 0
    return x, y


def balistic_simulated_bifurcation_gpu_naive(J, h, num_steps, time_step, num_trajectories: int, 
                                   a_0: Optional[float] = None, c_0_scaling: Optional[float] = None):
    if a_0 is None:
        a_0 = 1

    N, _ = J.shape
    mean_J = cp.sqrt(cp.sum(cp.square(J)) / (N * (N - 1)) )
    c_0 = 0.5 / (mean_J * sqrt(N))

    if c_0_scaling is not None:
        c_0 *= c_0_scaling


    a = cp.linspace(0, a_0, num=num_steps)

    x = cp.zeros((N, num_trajectories))
    y = cp.random.uniform(-0.1, 0.1, (N, num_trajectories))

    for t in tqdm(range(num_steps), desc="symulowana bifurkacja"):
        y += (-1 * (a_0 - a[t]) * x + c_0 * (cp.matmul(J,x) + h.reshape((N, 1)))) * time_step  # x(t)
        x += a_0 * y * time_step # x(t + 1)

        x, y = wall_gpu(x, y)

    x = cp.sign(x)
    return x, calculate_energy_gpu(J, h, x)

import matplotlib.pyplot as plt 
import numpy as np
from funkcje_pomocnicze import read_instance, test_pegasus, small_pegasus, P12

instance = P12

J, h = read_instance(instance.path, convention="minus_half")

n = len(h)
print(n)

J = cp.asarray(J, dtype=cp.float32)
h = cp.asarray(h, dtype=cp.float32)

states, energies = balistic_simulated_bifurcation_gpu_naive(J, h, time_step=0.5, num_steps=200, num_trajectories=2000, c_0_scaling=0.7)

# Konwersja z karty graficznej na CPU
energies = cp.asnumpy(energies)
states = cp.asnumpy(states)

print(f"Otrzymana energia: {min(energies)}")
print(f"Stan podstawowy: {instance.best_energy}")
print(f"Luka energetyczna: {((instance.best_energy - min(energies))/(2 * instance.best_energy) * 100):.2f}%")

unique_values, counts = np.unique(energies, return_counts=True)

2904

symulowana bifurkacja: 100%|██████████| 200/200 [00:13<00:00, 14.77it/s]

Otrzymana energia: -6713.0
Stan podstawowy: -6831.0
Luka energetyczna: 0.86%

# Implementacja SBM używając własnego kernelu

from setup import setup_paths
setup_paths()


import cupy as cp
from math import sqrt
from typing import Optional
from tqdm import tqdm

from funkcje_pomocnicze_gpu import calculate_energy_gpu

def balistic_simulated_bifurcation_gpu(J, h, num_steps, time_step, num_trajectories: int, 
                                   a_0: Optional[float] = None, c_0_scaling: Optional[float] = None):
    if a_0 is None:
        a_0 = 1
    
    N, _ = J.shape
    mean_J = np.sqrt(np.sum(np.square(J)) / (N * (N - 1)) )
    c_0 = 0.5 / (mean_J * sqrt(N))

    if c_0_scaling is not None:
        c_0 *= c_0_scaling

    dtype = cp.float32

    a = [dtype(i * a_0 / (num_steps - 1)) for i in range(num_steps)]
    
    a_0 = dtype(a_0)
    c_0 = dtype(c_0.item())
    time_step = dtype(time_step)

    x = cp.zeros((N, num_trajectories), dtype=dtype)
    y = cp.random.uniform(-0.1, 0.1, (N, num_trajectories), dtype=dtype)

    x_new = cp.empty_like(x)
    y_new = cp.empty_like(y)

    threadsperblock = 256  # Ilość wątków w bloku,
    blockspergrid_x = num_trajectories  # każdy blok zajmuje się trajektorią
    blockspergrid_y = (N + threadsperblock - 1) // threadsperblock  # wystarczająca ilość bloków by pomieścić całą kolumnę 
    blockspergrid = (blockspergrid_x, blockspergrid_y)
    
    kernel = cp.RawModule(path="cuda_kernels/sbm_kernel.ptx")
    update_y = kernel.get_function("update_y")
    update_x_and_wall = kernel.get_function("update_x_and_wall")

    for t in tqdm(range(num_steps), desc="balistyczna symulowana bifurkacja"):
        A = cp.matmul(J, x) + h.reshape((N, 1))
        update_y(blockspergrid, (threadsperblock,),(x, y, A, 
                                                    a_0, a[t], c_0, time_step, N,
                                                    y_new))
        y = y_new
        update_x_and_wall(blockspergrid, (threadsperblock,),(x, y, 
                                                             a_0, time_step, N, 
                                                             x_new, y_new))
        x = x_new
        y = y_new

    x = cp.sign(x)
    return x, calculate_energy_gpu(J, h, x)

# testy dla większej instancji

import matplotlib.pyplot as plt 
import numpy as np
from funkcje_pomocnicze import read_instance, test_pegasus, small_pegasus, P12

instance = P12

J, h = read_instance(instance.path, convention="minus_half")

J = cp.asarray(J, dtype=cp.float32)
h = cp.asarray(h, dtype=cp.float32)

states, energies = balistic_simulated_bifurcation_gpu(J, h, time_step=0.5, num_steps=200, num_trajectories=2000, c_0_scaling=0.7)

# Konwersja z karty graficznej na CPU
energies = cp.asnumpy(energies)
states = cp.asnumpy(states)

print(f"Otrzymana energia: {min(energies)}")
print(f"Stan podstawowy: {instance.best_energy}")
print(f"Luka energetyczna: {((instance.best_energy - min(energies))/(2 * instance.best_energy) * 100):.2f}%")

balistyczna symulowana bifurkacja: 100%|██████████| 200/200 [00:00<00:00, 6287.65it/s]

Otrzymana energia: -6713.0
Stan podstawowy: -6831.0
Luka energetyczna: 0.86%

# Dyskretna symulowana bifurkacja na gpu naiwna - zamieniamy np na cp

from setup import setup_paths
setup_paths()


import cupy as cp
from math import sqrt
from typing import Optional
from tqdm import tqdm

from funkcje_pomocnicze_gpu import calculate_energy_gpu


def wall_gpu(x: cp.ndarray, y: cp.ndarray):
    mask = cp.abs(x) > 1
    x[mask] = cp.sign(x[mask])
    y[mask] = 0
    return x, y


def discrete_simulated_bifurcation_gpu_naive(J, h, num_steps, time_step, num_trajectories: int, 
                                   a_0: Optional[float] = None, c_0_scaling: Optional[float] = None):
    if a_0 is None:
        a_0 = 1

    N, _ = J.shape
    mean_J = cp.sqrt(cp.sum(cp.square(J)) / (N * (N - 1)) )
    c_0 = 0.5 / (mean_J * sqrt(N))

    if c_0_scaling is not None:
        c_0 *= c_0_scaling


    a = cp.linspace(0, a_0, num=num_steps)

    x = cp.zeros((N, num_trajectories))
    y = cp.random.uniform(-0.1, 0.1, (N, num_trajectories))

    for t in tqdm(range(num_steps), desc="symulowana bifurkacja"):
        y += (-1 * (a_0 - a[t]) * x + c_0 * (cp.matmul(J,cp.sign(x)) + h.reshape((N, 1)))) * time_step  # x(t)
        x += a_0 * y * time_step # x(t + 1)

        x, y = wall_gpu(x, y)

    x = cp.sign(x)
    return x, calculate_energy_gpu(J, h, x)

# Dyskretna Symulowana Bifurkacja - własny kernel


from setup import setup_paths
setup_paths()

import cupy as cp
from math import sqrt
from typing import Optional
from tqdm import tqdm

from funkcje_pomocnicze_gpu import calculate_energy_gpu


def discrete_simulated_bifurcation_gpu(J, h, num_steps, time_step, num_trajectories: int, 
                                   a_0: Optional[float] = None, c_0_scaling: Optional[float] = None):
    if a_0 is None:
        a_0 = 1
    
    N, _ = J.shape
    mean_J = cp.sqrt(cp.sum(cp.square(J)) / (N * (N - 1)) )
    c_0 = 0.5 / (mean_J * sqrt(N))

    if c_0_scaling is not None:
        c_0 *= c_0_scaling

    dtype = cp.float32

    a = [dtype(i * a_0 / (num_steps - 1)) for i in range(num_steps)]
    
    a_0 = dtype(a_0)
    c_0 = dtype(c_0.item())
    time_step = dtype(time_step)

    x = cp.zeros((N, num_trajectories), dtype=dtype)
    y = cp.random.uniform(-0.1, 0.1, (N, num_trajectories), dtype=dtype)

    x_new = cp.empty_like(x)
    y_new = cp.empty_like(y)

    threadsperblock = 256  # Ilość wątków w bloku,
    blockspergrid_x = num_trajectories  # każdy blok zajmuje się trajektorią
    blockspergrid_y = (N + threadsperblock - 1) // threadsperblock  # wystarczająca ilość bloków by pomieścić całą kolumnę 
    blockspergrid = (blockspergrid_x, blockspergrid_y)
    
    kernel = cp.RawModule(path="cuda_kernels/sbm_kernel.ptx")
    update_y = kernel.get_function("update_y")
    update_x_and_wall = kernel.get_function("update_x_and_wall")

    for t in tqdm(range(num_steps), desc="dyskretna symulowana bifurkacja"):
        sign_x = cp.sign(x)
        A = cp.matmul(J, sign_x) + h.reshape((N, 1))
        update_y(blockspergrid, (threadsperblock,),(x, y, A, 
                                                    a_0, a[t], c_0, time_step, N,
                                                    y_new))
        y = y_new
        update_x_and_wall(blockspergrid, (threadsperblock,),(x, y, 
                                                             a_0, time_step, N,
                                                             x_new, y_new))
        x = x_new
        y = y_new

    x = cp.sign(x)
    return x, calculate_energy_gpu(J, h, x)

# testy dla większej instancji

import matplotlib.pyplot as plt 
import numpy as np
from funkcje_pomocnicze import read_instance, test_pegasus, small_pegasus, P8

instance = P8

J, h = read_instance(instance.path, convention="minus_half")

J = cp.asarray(J, dtype=cp.float32)
h = cp.asarray(h, dtype=cp.float32)

states, energies = discrete_simulated_bifurcation_gpu(J, h, time_step=0.5, num_steps=200, num_trajectories=2000, c_0_scaling=0.7)


print(f"Otrzymana energia: {min(energies)}")
print(f"Stan podstawowy: {instance.best_energy}")
print(f"Luka energetyczna: {((instance.best_energy - min(energies))/(2 * instance.best_energy) * 100):.2f}%")

dyskretna symulowana bifurkacja: 100%|██████████| 200/200 [00:00<00:00, 4180.57it/s]

Otrzymana energia: -2740.0
Stan podstawowy: -2752.0
Luka energetyczna: 0.22%

Wyżarzanie równoległe

# Implementacja algorytmu wyżarzania równoległego używając własnego kernelu w CUDA
from setup import setup_paths
setup_paths()


import os
import cupy as cp
from typing import Optional
from tqdm import tqdm

from funkcje_pomocnicze_gpu import calculate_energy_gpu



def parrarel_annealing_gpu(J, h, step_size: float, lambda_t_max: float, num_steps: int, num_trajectories: int,
                       schedule: Optional[list] = None, dtype = cp.float32):
    
    n = len(h)
    x = cp.zeros((n, num_trajectories), dtype=dtype)  # stan podstawowy dla H_innit = sum(x**2)
    momentum = cp.zeros((n, num_trajectories), dtype=dtype)
    state = cp.random.choice([dtype(-1.), dtype(1.)], size=(n, num_trajectories)) # losowy stan początkowy
    step_size = dtype(step_size)

    if schedule is None:
        schedule = [dtype(lambda_t_max * (1 - i/(num_steps -1))) for i in range(num_steps)]  # dlaczego nie linspace? chodzi o typy

    # cu_path = os.path.join("cuda_kernels", "pa_kernel.cu")
    # ptx_path = os.path.join("cuda_kernels", "pa_kernel.ptx")
    # command = f"nvcc --ptx --define-macro N={n} --define-macro M={num_trajectories} {cu_path} -o {ptx_path}"
    
    # os.system(command)
    kernel = cp.RawModule(path="cuda_kernels/pa_kernel.ptx")
    parrarel_annealing_step = kernel.get_function("parrarel_annealing_step")

    threadsperblock = 256  # Ilość wątków w bloku,
    blockspergrid_x = num_trajectories  # każdy blok zajmuje się trajektorią
    blockspergrid_y = (n + threadsperblock - 1) // threadsperblock  # wystarczająca ilość bloków by pomieścić całą kolumnę 
    blockspergrid = (blockspergrid_x, blockspergrid_y)

    x_new = cp.empty_like(x)
    momentum_new = cp.empty_like(momentum)
    state_new = cp.empty_like(state)

    for k in tqdm(range(num_steps), desc="wyżarzanie równoległe GPU"):
   
        lambda_t = schedule[k]
        A = cp.matmul(J, state)
        parrarel_annealing_step(blockspergrid, (threadsperblock,), (A, h, x, momentum, 
                                                               lambda_t, step_size, n,
                                                               momentum_new, x_new, state_new))
        momentum = momentum_new
        x = x_new
        state = state_new

    return state, calculate_energy_gpu(J, h, state)

from funkcje_pomocnicze import read_instance, small_pegasus, test_pegasus

instance = small_pegasus

J, h = read_instance(instance.path, convention="minus_half")


J = cp.asarray(J, dtype=cp.float32)
h = cp.asarray(h, dtype=cp.float32)

states, energies = parrarel_annealing_gpu(J, h, step_size=0.01, lambda_t_max=10, num_steps=1000, num_trajectories=5000)
print(f"Otrzymana energia: {min(energies)}")
print(f"Stan podstawowy: {instance.best_energy}")
print(f"Luka energetyczna: {((instance.best_energy - min(energies))/(2 * instance.best_energy) * 100):.2f}%")

wyżarzanie równoległe GPU: 100%|██████████| 1000/1000 [00:00<00:00, 9518.56it/s]

Otrzymana energia: -39.0
Stan podstawowy: -39.0
Luka energetyczna: -0.00%

Bruteforce

Ta implementacja pochodzi z artykułu:

Jałowiecki, K., Rams, M. M., & Gardas, B. (2021). Brute-forcing spin-glass problems with CUDA. Computer Physics Communications, 260, 107728. DOI: 10.1016/j.cpc.2020.107728

# brute force CUDA

import cupy as cp
from tqdm import tqdm

def select_lowest(energies, states, num_states):
    
    indices = cp.argpartition(energies, num_states)[:num_states]

    low_energies = energies[indices]
    low_states = states[indices]
    return low_energies, low_states

def sort_by_key(energies, states, num_states):
    order = cp.argsort(energies)[:num_states]

    low_energies = energies[order]
    low_states = states[order]
    return low_energies, low_states


def brute_force_gpu(Q,  num_states: int, sweep_size_exponent: int = 10, threadsperblock: int = 256):
    N, _ = Q.shape
    if N > 64:
        raise ValueError("Za wysoka wartość N. Ta implementacja wspiera co najwyżej 64 spiny (64-bitowy integer)")
    sweep_size = 2**sweep_size_exponent
    num_chunks = 2**(N-sweep_size_exponent)

   
    brute_force_kernel = cp.RawModule(path="cuda_kernels/brute_force_kernel.ptx")
    compute_energies = brute_force_kernel.get_function("compute_energies")

    blockspergrid = sweep_size//threadsperblock
    
    final_energies = cp.array([])
    final_states = cp.array([])
    
    for i in tqdm(range(num_chunks), desc="wyczerpujące przeszukiwanie"):

        energies = cp.empty(sweep_size, dtype=cp.float32)
        states = cp.empty(sweep_size, dtype=cp.int64)
        compute_energies((blockspergrid,), (threadsperblock,), 
                         (Q, cp.int32(N), cp.int32(sweep_size_exponent), energies, states, cp.int64(i)))

        low_energies, low_states = select_lowest(energies, states, num_states)
        
        final_energies = cp.concatenate((final_energies, low_energies))
        final_states = cp.concatenate((final_states, low_states))
        if i != 0:

            final_energies, final_states = select_lowest(final_energies, final_states, num_states)

    return sort_by_key(final_energies, final_states, num_states)

    

# test
import time
import os
import cupy as cp

from tqdm import tqdm
from IPython.utils.io import capture_output

n_range = list(range(15, 32))

# Kompilacja
Q = cp.random.uniform(-1, 1, size=(10, 10), dtype=cp.float32)
Q = cp.triu(Q)
with capture_output() as captured:
    brute_force_gpu(Q, 10)

times = []
for n in tqdm(n_range, desc="Wyczerpujące przeszukiwanie dla różnych n"):
    Q = cp.random.uniform(-1, 1, size=(n, n), dtype=cp.float32)
    Q = cp.triu(Q)
    start = time.time()
    with capture_output() as captured:
        brute_force_gpu(Q, num_states=10, sweep_size_exponent=min(n, 21), threadsperblock=512)
    end = time.time()
    times.append(end - start)



plt.figure(figsize=(8, 5))
plt.plot(n_range, times, marker='o', linestyle='-', color='blue', label='Execution Time')
plt.xlabel('Rozmiar (ilość spinów)')
plt.xticks(n_range)
plt.yscale('log', base=10)
plt.ylabel('Czas obliczeń (s)')
plt.grid(True)
plt.show()

Wyczerpujące przeszukiwanie dla różnych n: 100%|██████████| 17/17 [00:03<00:00,  4.56it/s]