如何有效计算核/矩阵

Question

import numpy as np

solution_point_count = 150
impedance_datapoint_count = 134 

impedance_frequency = np.logspace(np.log10(100000), np.log10(0.0199), impedance_datapoint_count)
solution_frequency_logspace = np.logspace(np.log10(100000), np.log10(0.06), solution_point_count)
    
# Initializing zeros matrix to store results     
Z_RC = np.zeros((impedance_datapoint_count,solution_point_count))

#Constructing Debye Model for RC Kernel
#Debye Model: Z_RC = 1/(1 + iω*𝜏) where ω = 2πf and 𝜏 = RC = 1/(2πf)

for n in range(solution_point_count):
    for m in range(impedance_datapoint_count):
        Equation = 1/(1+1j*impedance frequency[m]/solution_frequency_logspace[n])
            Z_RC[m,n] = np.real(Equation) # Extracting real part only

虽然这可以创建内核 Z_RC，但它非常慢并且可能效率低下。有没有更有效的方法来做到这一点？也许通过矢量化？

Answer 1

最简单的解决方案是使用 numba:

from numba import njit

@njit
def calculate_numba(impedance_frequency, solution_frequency_logspace):
    solution_point_count = len(solution_frequency_logspace)
    impedance_datapoint_count = len(impedance_frequency)

    Z_RC = np.empty((impedance_datapoint_count, solution_point_count))

    for n in range(solution_point_count):
        for m in range(impedance_datapoint_count):
            Equation = 1 / (
                1 + 1j * impedance_frequency[m] / solution_frequency_logspace[n]
            )
            Z_RC[m, n] = np.real(Equation)  # Extracting real part only

    return Z_RC

---

使用

perfplot

进行基准测试：

import numpy as np
import perfplot
from numba import njit


def generate_impedance_frequency(solution_point_count, impedance_datapoint_count):
    impedance_frequency = np.logspace(
        np.log10(100000), np.log10(0.0199), impedance_datapoint_count
    )
    solution_frequency_logspace = np.logspace(
        np.log10(100000), np.log10(0.06), solution_point_count
    )
    return impedance_frequency, solution_frequency_logspace


def calculate_normal(impedance_frequency, solution_frequency_logspace):
    solution_point_count = len(solution_frequency_logspace)
    impedance_datapoint_count = len(impedance_frequency)

    Z_RC = np.empty((impedance_datapoint_count, solution_point_count))

    for n in range(solution_point_count):
        for m in range(impedance_datapoint_count):
            Equation = 1 / (
                1 + 1j * impedance_frequency[m] / solution_frequency_logspace[n]
            )
            Z_RC[m, n] = np.real(Equation)  # Extracting real part only

    return Z_RC


@njit
def calculate_numba(impedance_frequency, solution_frequency_logspace):
    solution_point_count = len(solution_frequency_logspace)
    impedance_datapoint_count = len(impedance_frequency)

    Z_RC = np.empty((impedance_datapoint_count, solution_point_count))

    for n in range(solution_point_count):
        for m in range(impedance_datapoint_count):
            Equation = 1 / (
                1 + 1j * impedance_frequency[m] / solution_frequency_logspace[n]
            )
            Z_RC[m, n] = np.real(Equation)  # Extracting real part only

    return Z_RC


# check if the calculation is the same:

solution_point_count = 150
impedance_datapoint_count = 134

impedance_frequency, solution_frequency_logspace = generate_impedance_frequency(
    solution_point_count, impedance_datapoint_count
)

out1 = calculate_normal(impedance_frequency, solution_frequency_logspace)
out2 = calculate_numba(impedance_frequency, solution_frequency_logspace)

assert np.allclose(out1, out2)


perfplot.show(
    setup=lambda n: generate_impedance_frequency(n, n),
    kernels=[
        calculate_normal,
        calculate_numba,
    ],
    labels=["normal", "numba"],
    n_range=[10, 20, 50, 100, 150, 200],
    xlabel="N x N",
    logx=True,
    logy=True,
)

结果：