如何在Python中并行运行SMO算法中的两个核矩阵计算任务？

Great question! Let's break this down clearly: you can't directly parallelize Kik = self.kernel_matrix(X, i) and Kjk = self.kernel_matrix(X, j) in your current code structure—because Kjk depends on finding j first, and j is computed after Kik in your existing flow.

But here's the good news: the computation of Kik and the process of finding j are completely independent of each other. We can parallelize those two tasks instead, which will still give you a meaningful speedup, especially if your kernel calculations are computationally heavy (like RBF kernels on large datasets).

Why This Works #

In your current code, the sequence is:

1. Find index i
2. Compute Kik (slow kernel matrix operation)
3. Calculate index j (fast array operations)
4. Compute Kjk (another slow kernel matrix operation)

Since step 2 (computing Kik) doesn't rely on anything from step 3 (finding j), and vice versa, we can run them at the same time. Once both finish, we can compute Kjk normally.

How to Modify the Code #

We'll use Python's concurrent.futures.ThreadPoolExecutor to run the independent tasks in parallel. Here's the step-by-step adjustment:

Step 1: Import the Required Module #

Add this at the top of your script:

from concurrent.futures import ThreadPoolExecutor

Step 2: Refactor for Parallel Execution #

Replace this sequential block in your smo method:

i = numpy.argmax(yg_i)
Kik = self.kernel_matrix(X, i)
indices_violating_Ai_1 = y_pos_ind * alpha_neg_ind
indices_violating_Ai_2 = y_neg_ind * alpha_pos_ind
indices_violating_Ai = indices_violating_Ai_1 + indices_violating_Ai_2
yg_j = yg.copy()
yg_j[indices_violating_Ai] = float('+inf')
# Second of the maximum violating pair
j = numpy.argmin(yg_j)
Kjk = self.kernel_matrix(X, j)

With this parallelized version:

i = numpy.argmax(yg_i)

# Helper function to compute j (needs to be callable for the executor)
def find_j(y_pos_ind, alpha_neg_ind, y_neg_ind, alpha_pos_ind, yg):
    indices_violating_Ai_1 = y_pos_ind * alpha_neg_ind
    indices_violating_Ai_2 = y_neg_ind * alpha_pos_ind
    indices_violating_Ai = indices_violating_Ai_1 + indices_violating_Ai_2
    yg_j = yg.copy()
    yg_j[indices_violating_Ai] = float('+inf')
    return numpy.argmin(yg_j)

# Run Kik computation and j-finding in parallel
with ThreadPoolExecutor(max_workers=2) as executor:
    future_kik = executor.submit(self.kernel_matrix, X, i)
    future_j = executor.submit(find_j, y_pos_ind, alpha_neg_ind, y_neg_ind, alpha_pos_ind, yg)
    
    # Retrieve results once both tasks complete
    Kik = future_kik.result()
    j = future_j.result()

# Compute Kjk now that we have j
Kjk = self.kernel_matrix(X, j)

Key Considerations #

• ThreadPool vs ProcessPool: If your kernel_matrix uses pure Python code (not optimized numpy operations that release the GIL), switch to ProcessPoolExecutor instead. Just ensure all arguments passed to the executor are serializable (numpy arrays and basic types work fine).
• Overhead Check: Parallelization adds a tiny amount of overhead. This is only worth it if your kernel_matrix calls are slow. For small datasets or linear kernels, the overhead might cancel out any speed gain.
• Algorithm Correctness: This refactor doesn't change the SMO algorithm's logic—you'll get the exact same alpha results as the sequential version, just potentially faster.

Full Modified smo Method #

import numpy
from concurrent.futures import ThreadPoolExecutor

def smo(self, X, y):
    iterations = 0
    n_samples = X.shape[0]
    # Initial coefficients
    alpha = numpy.zeros(n_samples)
    # Initial gradient
    g = numpy.ones(n_samples)

    while True:
        yg = g * y
        # KKT Conditions
        y_pos_ind = (y == 1)
        y_neg_ind = (numpy.ones(n_samples) - y_pos_ind).astype(bool)
        alpha_pos_ind = (alpha >= self.C)
        alpha_neg_ind = (alpha <= 0)

        indices_violating_Bi_1 = y_pos_ind * alpha_pos_ind
        indices_violating_Bi_2 = y_neg_ind * alpha_neg_ind
        indices_violating_Bi = indices_violating_Bi_1 + indices_violating_Bi_2

        yg_i = yg.copy()
        yg_i[indices_violating_Bi] = float('-inf')
        # First of the maximum violating pair
        i = numpy.argmax(yg_i)

        # --- Parallelized Block Start ---
        def find_j(y_pos_ind, alpha_neg_ind, y_neg_ind, alpha_pos_ind, yg):
            indices_violating_Ai_1 = y_pos_ind * alpha_neg_ind
            indices_violating_Ai_2 = y_neg_ind * alpha_pos_ind
            indices_violating_Ai = indices_violating_Ai_1 + indices_violating_Ai_2
            yg_j = yg.copy()
            yg_j[indices_violating_Ai] = float('+inf')
            return numpy.argmin(yg_j)

        with ThreadPoolExecutor(max_workers=2) as executor:
            future_kik = executor.submit(self.kernel_matrix, X, i)
            future_j = executor.submit(find_j, y_pos_ind, alpha_neg_ind, y_neg_ind, alpha_pos_ind, yg)
            Kik = future_kik.result()
            j = future_j.result()
        # --- Parallelized Block End ---

        Kjk = self.kernel_matrix(X, j)

        # Optimality criterion
        if(yg_i[i] - yg_j[j]) < self.tol or (iterations >= self.max_iter):
            break

        min_term_1 = (y[i] == 1) * self.C - y[i] * alpha[i]
        min_term_2 = y[j] * alpha[j] + (y[j] == -1) * self.C
        min_term_3 = (yg_i[i] - yg_j[j]) / (Kik[i] + Kjk[j] - 2 * Kik[j])
        # Direction search
        lamda = numpy.min([min_term_1, min_term_2, min_term_3])

        # Gradient update
        g += lamda * y * (Kjk - Kik)
        # Update coefficients
        alpha[i] = alpha[i] + y[i] * lamda
        alpha[j] = alpha[j] - y[j] * lamda

        iterations += 1

    print('{} iterations to arrive at the minimum'.format(iterations))
    return alpha

内容的提问来源于stack exchange，提问作者Shankar kumar