我想实现一些快速的凸分析运算--近似算子等等.我是Cython的新手,我认为这将是适合这项工作的工具.我有纯Python和Cython语言(下面是mwe_py.pymwe_c.pyx)的实现.然而,当我比较它们时,Python+Numpy版本比Cython版本要快得多.这是为什么?我试过使用内存视图,这应该允许更快的索引/操作;然而,性能差异非常明显!任何关于如何修复下面的mwe_c.pyx以接近"最佳"Cython代码的建议都将不胜感激.

import pyximport; pyximport.install(language_level=3)

import mwe_c
import mwe_py
import numpy as np
from time import time

n = 100000
nreps = 10000
x = np.random.randn(n)
z = np.random.randn(n)
tau = 1.0

t0 = time()
for _ in range(nreps):
    out = mwe_c.prox_translation(mwe_c.prox_two_norm, x, z, tau)
t1 = time()
print(t1 - t0)


t0 = time()
for _ in range(nreps):
    out = mwe_py.prox_translation(mwe_py.prox_two_norm, x, z, tau)
t1 = time()
print(t1 - t0)

它分别给出了以下输出:

10.76103401184082  # (seconds)
5.988733291625977  # (seconds)

下面是mwe_py.py:

import numpy.linalg as la

def proj_two_norm(x):
    """projection onto l2 unit ball"""
    X = la.norm(x)
    if X <= 1:
        return x
    return x / X

def prox_two_norm(x, tau):
    """proximal mapping of l2 norm with parameter tau"""
    return x - proj_two_norm(x / tau)

def prox_translation(prox_func, x, z, tau=None):
  """Compute prox(f(. - z))(x) where prox_func(x, tau) is prox(tau * f)(x)."""
  if tau is None:
      tau = 1.0
  return z + prox_func(x - z, tau)

下面是最后的mwe_c.pyx条:

import numpy as np
cimport numpy as np


cdef double [::1] aasubtract(double [::1] x, double [::1] z):
    cdef unsigned int i, m = len(x), n = len(z);
    assert m == n, f"vectors must have the same length"
    cdef double [::1] out = np.copy(x);
    for i in range(n):
        out[i] -= z[i]
    return out


cdef double [::1] vsdivide(double [::1] x, double tau):
    """Divide an array by a scalar element-wise."""
    cdef:
        unsigned int i, n = len(x);
        double [::1] out = np.copy(x);
    for i in range(n):
        out[i] /= tau
    return out


cdef double two_norm(double [::1] x):
    cdef:
        double out = 0.0;
        unsigned int i, n=len(x);
    for i in range(n):
        out = out + x[i]**2
    out = out **.5
    return out


cdef double [::1] proj_two_norm(double [::1] x):
    """project x onto the unit two ball."""
    cdef double x_norm = two_norm(x);
    cdef unsigned int i, n = len(x);
    cdef double [::1] p = np.copy(x);
    if x_norm <= 1:
        return p
    for i in range(n):
        p[i] = p[i] / x_norm
    return p


cpdef double [::1] prox_two_norm(double [::1] x, double tau):
    """double [::1] prox_two_norm(double [::1] x, double tau)"""
    cdef unsigned int i, n = len(x);
    cdef double [::1] out = x.copy(), Px = x.copy();
    Px = proj_two_norm(vsdivide(Px, tau));
    for i in range(n):
        out[i] = out[i] - Px[i]
    return out


cpdef prox_translation(
    prox_func,
    double [::1] x,
    double [::1] z,
    double tau=1.0
):
    cdef:
        unsigned int i, n = len(x);
        double [::1] out = prox_func(aasubtract(x, z), tau);
    for i in range(n):
        out[i] += z[i];
    return out

推荐答案

主要问题是,您将优化后的Numpy代码与优化程度较低的Cython代码进行比较.事实上,Numpy利用了SIMD instructions(就像x86-64处理器上的SSE和AVX/AVX2),它们能够连续计算许多项.默认情况下,Cython使用默认的-O2优化级别,这不会启用任何自动向量化策略,从而导致标量代码变慢(除非您使用最新版本的GCC).您可以使用-O3来告诉大多数编译器(例如,老GCC和Cang),以启用自动矢量化.请注意,这不足以产生非常快的代码.事实上,出于兼容性的考虑,编译器只在x86-64处理器上使用传统的SIMD指令.-mavx-mavx2启用AVX/AVX-2指令集,以便生成更快的代码,假设您的机器支持它(否则它只会崩溃).-mfma可能也会有所帮助.-march=native也可用于 Select 目标平台上可用的最佳指令集.请注意,Numpy在运行时(部分地)执行此判断(这要归功于GCC特定的C特性).

第二个主要问题是,如果不违反IEEE-754标准,则out = out + x[i]**2会产生一个很大的loop-carried dependency chain,编译器无法对其进行优化.事实上,有很长的加法链要执行,处理器执行加法的速度不可能比使用当前代码串行执行每条加法指令更快.问题是,将两个浮点数相加会有相当大的延迟(在相当现代的x86-64处理器上通常是3到4个周期).这意味着处理器不能pipeline该指令.事实上,现代处理器通常可以(每个核)并行执行两个加法,但电流环路阻止了这一点.最后,这个循环完全是latency-bound.您可以手动将此问题修复为unrolling the loop.

Using -ffast-math can help compilers to do such optimizations but at the expense of breaking the IEEE-754 standard. If you use this option, then you need not to use special values like NaN numbers nor some operations. For more information, please read What does gcc's ffast-math actually do? .

此外,请注意array copies are expensive份,我不确定是否需要所有的副本.您可以创建一个新的空数组并填充它,而不是对数组的副本进行操作.这将会更快,特别是对于大型array.

最后,divisions are slow条.请考虑乘以倒数.由于IEEE-754标准,编译器无法执行此优化,但您可以轻松地执行此操作.也就是说,您需要确保这在您的情况下是可以的,因为它可以稍微改变结果.使用-ffast-math也应该会自动修复此问题.

请注意,Numpy许多开发人员都知道编译器和处理器是如何工作的,所以他们已经进行了手动优化,以便生成快速的代码(就像我多次做的那样).在处理大型数组时,除非您合并循环或使用多线程,否则很难胜过Numpy.事实上,与计算单元和Numpy Create My temporary arrays相比,如今的RAM相当慢.可以使用Cython来避免创建大多数这样的临时array.

Python相关问答推荐

如何从pandas DataFrame中获取. groupby()和. agg()之后的子列?

使用__json__的 pyramid 在客户端返回意外格式

pandas fill和bfill基于另一列中的条件

jsonschema日期格式

如何防止html代码出现在quarto gfm报告中的pandas表之上

Django.core.exceptions.SynchronousOnlyOperation您不能从异步上下文中调用它-请使用线程或SYNC_TO_ASYNC

在matplotlib中重叠极 map 以创建径向龙卷风图

具有不同坐标的tkinter canvs.cocords()和canvs.moveto()

我的浮点问题--在C++/Python中的试用

捕获脚本和退出代码的多行输出

基于符号和位置 Select 数据帧特定区域的毕达式方法

django中没有预期的输出

PANDA TO_DICT-按键列出行(_D)

允许在枚举中使用一组特定的未定义值

是否 Select 所有整型列,除了几个python-polars列?

支持向量机模型突出错误的数据点作为支持向量

POLARS:从GROUP_BY列表中 Select 值,并从另一列中 Select 值

将Hangman游戏中的&替换为所有比赛的玩家猜测

Pandas 多列数据帧的重采样和内插

使用Python下载文件时出现错误699