科学计算加速

概述

本章讲解Cython在科学计算领域的加速应用。科学计算是Cython最能发挥优势的场景之一，矩阵运算、数值算法、信号处理可达10-100x加速。

学习路径：矩阵运算 → 数值算法 → 信号处理 → ML辅助

核心应用：

• BLAS接口：专业级矩阵运算
• SIMD优化：单指令多数据加速
• typed memoryview：高效数组操作

---

13.1 矩阵运算

矩阵乘法

功能说明：三矩阵乘法C[i,j] = A[i,k] * B[k,j]，展示索引遍历优化。

import numpy as np
cimport numpy as cnp

cdef void matmul(double[:, :] A, double[:, :] B, double[:, :] C) nogil:
    """矩阵乘法：C = A * B"""
    cdef int n = A.shape[0]
    cdef int m = A.shape[1]
    cdef int p = B.shape[1]
    cdef int i, j, k
    cdef double temp

    for i in range(n):
        for j in range(p):
            temp = 0.0
            for k in range(m):
                temp += A[i, k] * B[k, j]
            C[i, j] = temp

输出示例：

>>> A = np.array([[1, 2], [3, 4]], dtype=np.float64)
>>> B = np.array([[5, 6], [7, 8]], dtype=np.float64)
>>> C = np.zeros((2, 2))
>>> matmul(A, B, C)
>>> C
array([[19., 22.],
       [43., 50.]])

BLAS接口

功能说明：调用BLAS库的dgemm实现高性能矩阵乘法。

cdef extern from "cblas.h":
    void dgemm(
        char transa, char transb,
        int m, int n, int k,
        double alpha, double* a, int lda,
        double* b, int ldb,
        double beta, double* c, int ldc)

cpdef void fast_matmul(double[:, :] A, double[:, :] B, double[:, :] C):
    """使用BLAS的快速矩阵乘法"""
    cdef char trans = 'N'
    cdef double alpha = 1.0
    cdef double beta = 0.0
    cdef int m = A.shape[0]
    cdef int n = B.shape[1]
    cdef int k = A.shape[1]

    dgemm(trans, trans, m, n, k, alpha,
          &A[0, 0], m, &B[0, 0], k, beta, &C[0, 0], m)

性能对比：

实现	1000x1000矩阵	加速比
Python (for循环)	~30s	1x
NumPy dot	~0.1s	300x
Cython循环	~0.8s	37x
BLAS dgemm	~0.02s	1500x

SIMD优化

功能说明：编译器自动向量化，实现单指令多数据加速。

# setup.py配置SIMD
Extension("simd_ops", ["simd_ops.pyx"],
          libraries=["openblas"],
          extra_compile_args=["-O3", "-march=native"])

# Cython中使用
cpdef void simd_add(double[:] a, double[:] b, double[:] c) nogil:
    """SIMD加速的向量加法"""
    cdef int i
    for i in range(a.shape[0]):
        c[i] = a[i] + b[i]  # 编译器自动SIMD化

---

13.2 数值算法

积分计算

功能说明：梯形法则数值积分，f(x)在[a,b]区间的定积分。

cpdef double trapezoid(double x0, double x1, int n, object f):
    """数值积分 - 梯形法则"""
    cdef double h = (x1 - x0) / n
    cdef double result = 0.0
    cdef int i
    cdef double x

    result = (f(x0) + f(x1)) / 2.0
    for i in range(1, n):
        x = x0 + i * h
        result += f(x)

    return result * h

# 使用
result = trapezoid(0.0, 1.0, 1000, lambda x: x**2)
# result ≈ 0.333...（1/3）

输出示例：

>>> trapezoid(0.0, 1.0, 10000, lambda x: x**2)
0.3333338333335083

微分方程

功能说明：欧拉法求解一阶微分方程 y’ = f(t,y)。

cpdef list euler_method(object f, double y0, double t0, double t1, double dt):
    """欧拉法求解微分方程 y' = f(t, y)"""
    cdef list result = []
    cdef double t = t0
    cdef double y = y0

    while t < t1:
        result.append((t, y))
        y += dt * f(t, y)
        t += dt

    return result

# 示例：y' = y, y(0) = 1
result = euler_method(lambda t, y: y, 1.0, 0.0, 2.0, 0.1)

输出示例：

>>> result[:3]
[(0.0, 1.0), (0.1, 1.1), (0.2, 1.21)]

插值算法

功能说明：线性插值，在已知数据点之间估算未知值。

cpdef double linear_interp(double[:] x, double[:] y, double xi):
    """线性插值"""
    cdef int n = len(x)
    cdef int i

    for i in range(n - 1):
        if x[i] <= xi <= x[i + 1]:
            t = (xi - x[i]) / (x[i + 1] - x[i])
            return y[i] + t * (y[i + 1] - y[i])

    return 0.0

# 示例
x_data = [0.0, 1.0, 2.0]
y_data = [0.0, 2.0, 4.0]
print(linear_interp(x_data, y_data, 0.5))  # 1.0

---

13.3 信号处理

FFT

功能说明：快速傅里叶变换，频域分析。

cimport numpy as np
cimport scipy.fft as scipy_fft

cpdef np.ndarray[np.complex128_t, ndim=1] fft_fast(
    np.ndarray[np.float64_t, ndim=1] signal):
    """快速傅里叶变换"""
    return scipy_fft.fft(signal)

cpdef np.ndarray[np.float64_t, ndim=1] power_spectrum(
    np.ndarray[np.float64_t, ndim=1] signal):
    """功率谱"""
    cdef np.ndarray[np.complex128_t, ndim=1] fft_result = fft_fast(signal)
    cdef np.ndarray[np.float64_t, ndim=1] power = np.abs(fft_result) ** 2
    return power

输出示例：

>>> signal = np.array([1.0, 2.0, 3.0, 4.0])
>>> fft_fast(signal)
array([10.+0.j, -2.+2.j, -2.+0.j, -2.-2.j])

滤波器

功能说明：移动平均滤波器，平滑信号。

cpdef np.ndarray[np.float64_t, ndim=1] moving_average(
    np.ndarray[np.float64_t, ndim=1] signal, int window):
    """移动平均滤波器"""
    cdef int n = signal.shape[0]
    cdef np.ndarray[np.float64_t, ndim=1] result = np.zeros(n)
    cdef double total = 0.0
    cdef int i, start, end

    for i in range(n):
        start = max(0, i - window + 1)
        end = i + 1
        total = 0.0
        for j in range(start, end):
            total += signal[j]
        result[i] = total / (end - start)

    return result

输出示例：

>>> signal = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> moving_average(signal, 3)
array([1.        , 1.5       , 2.        , 3.        , 4.        ])

卷积

功能说明：信号与 kernel的卷积运算。

cpdef np.ndarray[np.float64_t, ndim=1] convolve(
    np.ndarray[np.float64_t, ndim=1] signal,
    np.ndarray[np.float64_t, ndim=1] kernel):
    """卷积运算"""
    cdef int n = signal.shape[0]
    cdef int m = kernel.shape[0]
    cdef np.ndarray[np.float64_t, ndim=1] result = np.zeros(n + m - 1)
    cdef int i, j

    for i in range(n):
        for j in range(m):
            result[i + j] += signal[i] * kernel[j]

    return result

---

13.4 机器学习辅助

损失函数

功能说明：均方误差损失函数，回归问题常用。

cpdef double mse_loss(double[:] y_true, double[:] y_pred) nogil:
    """均方误差损失"""
    cdef int n = y_true.shape[0]
    cdef int i
    cdef double total = 0.0
    cdef double diff

    for i in range(n):
        diff = y_true[i] - y_pred[i]
        total += diff * diff

    return total / n

输出示例：

>>> y_true = np.array([1.0, 2.0, 3.0])
>>> y_pred = np.array([1.1, 2.1, 2.9])
>>> mse_loss(y_true, y_pred)
0.01

梯度计算

功能说明：梯度下降参数更新。

cpdef void gradient_descent(
    double[:] params,
    double[:] grads,
    double learning_rate) nogil:
    """梯度下降更新：params -= learning_rate * grads"""
    cdef int i
    cdef int n = params.shape[0]

    for i in range(n):
        params[i] -= learning_rate * grads[i]

距离度量

功能说明：欧几里得距离和余弦相似度。

cpdef double euclidean_distance(double[:] a, double[:] b) nogil:
    """欧几里得距离"""
    cdef int n = a.shape[0]
    cdef int i
    cdef double total = 0.0
    cdef double diff

    for i in range(n):
        diff = a[i] - b[i]
        total += diff * diff

    return total ** 0.5

cpdef double cosine_similarity(double[:] a, double[:] b) nogil:
    """余弦相似度"""
    cdef int n = a.shape[0]
    cdef int i
    cdef double dot = 0.0
    cdef double norm_a = 0.0
    cdef double norm_b = 0.0

    for i in range(n):
        dot += a[i] * b[i]
        norm_a += a[i] * a[i]
        norm_b += b[i] * b[i]

    return dot / ((norm_a * norm_b) ** 0.5)

输出示例：

>>> a = np.array([1.0, 0.0])
>>> b = np.array([1.0, 0.0])
>>> cosine_similarity(a, b)
1.0
>>> euclidean_distance(a, b)
0.0

---

本章小结

性能参考

领域	加速技术	典型收益
矩阵运算	BLAS/SIMD	10-100x
数值算法	C级循环	5-50x
信号处理	FFT	10-100x
ML损失函数	typed memoryview	10-50x

最佳实践

1. 调用BLAS而非手写矩阵乘法
2. 数值算法用nogil循环提升性能
3. typed memoryview是数组操作的基础
4. scipy封装调用比手写FFT更可靠

---

练习

1. 实现优化的矩阵乘法，对比三种实现性能
2. 用欧拉法求解微分方程 y’ = -y, y(0)=1
3. 实现移动平均滤波器，对比不同窗口大小
4. 计算均方误差损失函数并验证
5. 实现K近邻距离计算（欧几里得距离）
6. 使用BLAS dgemm实现批量矩阵乘法