Part1. 加载/存储组的向量化

GCC 5.0 显著的提升了 vector
向量的加载和存储组的代码质量,我这里说的是连续顺序的迭代,例如:

x = a[i], y = a[i + 1], z = a[i + 2] 通过 i 进行迭代,加载了大小为
3 的组

组大小由加载和存储地址的最大和最小值来确定,例如 (i + 2) – (i) + 1 = 3

组中加载和存储的次数小于和等于组的大小,例如:

x = a[i], z = a[I + 2] 通过 i 进行迭代,尽管只有 2
次加载,但是加载组的大小为 3

GCC 4.9 向量组的大小是 2 的指数,而 GCC 5.0 向量化组的大小是 3
,也可以是 2 的指数,其他的组大小使用比较少。

最常用加载和存储组的场景是结构数组。

  1. 图像转换 (例如将 RGB 结构转为其他)

    (场景测试 )

  2. 多维的坐标 (测试场景
    )

  3. 向量的乘法常数矩阵

a[i][0] = 7 * b[i][0] – 3 * b[i][1];
a[i][1] = 2 * b[i][0] + b[i][1];

基本上 GCC 5.0 给我们带来了:

  1. 引入大小为 3 的向量加载和存储组

  2. 改进对原有支持的其他组大小

  3. 通过为特定的 x86 CPU 优化的代码来最大化加载和存储组性能

下面是一个用来比较 GCC 4.9 和 GCC 5.0
性能的一段代码(最大化向量中的元素个数)

?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int i, j, k; 
byte *in = a, *out = b;
for (i = 0; i < 1024; i++)
{
  for (k = 0; k < STGSIZE; k++)
    {
      byte s = 0;
      for (j = 0; j < LDGSIZE; j++)
        s += in[j] * c[j][k];
      out[k] = s;
    }
  in += LDGSIZE;
  out += STGSIZE;
}

而 “c” 是一个固定的矩阵:

?

1
2
3
4
5
6
7
8
const byte c[8][8] = {1, -1, 1, -1, 1, -1, 1, -1,
                      1, 1, -1, -1, 1, 1, -1, -1,
                      1, 1, 1, 1, -1, -1, -1, -1,
                      -1, 1, -1, 1, -1, 1, -1, 1,
                      -1, -1, 1, 1, -1, -1, 1, 1,
                      -1, -1, -1, -1, 1, 1, 1, 1,
                      -1, -1, -1, 1, 1, 1, -1, 1,
                      1, -1, 1, 1, 1, -1, -1, -1};

在循环中使用简单的计算,例如 add、sub 等操作速度非常快。

  • “in“ 和 “out” 指针指向全局数组 “a[1024 * LDGSIZE]” 和 “b[1024 *
    STGSIZE]”

  • byte 是一个无符号的 char

  • LDGSIZE 和 STGSIZE – 根据组大小定义加载和存储组的宏

编译选项 “-Ofast” 加 “-march=slm” 用于 Silvermont, “-march=core-avx2”
用于 Haswell 以及所有合并 -DLDGSIZE={1,2,3,4,8} -DSTGSIZE={1,2,3,4,8}

GCC 5.0 到 4.9 (时间上,越大越好)

Silvermont: Intel(R) Atom(TM) CPU  C2750  @ 2.41GHz

性能提升 6.5 倍

图片 1

上 表结果我们可以看到组大小为 3 的时候结果没那么好。因为当组大小为 3
时在 Slivermont 上需要 8 个 pshufb 指令和大约 5 个
ticks。当然循环依然是向量化的,如果在循环内有更消耗 CPU
的计算,那么效果还是会很不错。(我们再看其他组大小)

Haswell: Intel(R) Core(TM) i7-4770K CPU @ 3.50GHz

性能提升 3 倍!

图片 2

在 Haswell 上组大小为 3 的时候只需 1 个
tick。我们能看到最大的提升就是组大小为 3 的时候。

上述的实验你可以通过下面地址获取相应编译器

GCC 4.9:

GCC 5.0 trunk built at revision 217914:

Download
图片 3matrix.zip

via
intel

(文/开源中国)    

可以自动矢量化的代码


使用gcc version 4.3.4的gcc编译器,编译下面的示例代码:

int main(){

  const unsigned int ArraySize = 10000000;

  float* a = new float[ArraySize];

  float* b = new float[ArraySize];

  float* c = new float[ArraySize];

  for (unsigned int j = 0; j< 200 ; j++) // some repetitions

    for ( unsigned int i = 0; i < ArraySize; ++ i)

        c[i] = a[i] * b[i];

}

编译并运行上面的代码:

g++ vectorisable_example.cpp -o vectorisable_example -O2 time ./vectorisable_example

简单,现在我们看下如何告诉 compiler 自动对代码进行矢量化:

g++ vectorisable_example.cpp -o vectorisable_example_vect -O2 -ftree-vectorize time ./vectorisable_example_vect

第二种编译方式,产生的可执行程序比第一种产生的可执行程序运行的更快?我们深入研究下,看看为什么他就这么快?

1 向量化( Vectorization )

简介


从上世纪九十年代后期开始,英特尔就已经将单指令多数据(SIMD)机器指令集成到其商品CPU产品线中。这些SIMD指令允许一次将相同的数学运算(乘,加…)应用于2,4或甚至更多的数据值上。这组数据值也被称为“vector”(不要与代数中的vector混淆)。理论上矢量计算的性能增益等于CPU可以容纳的矢量单位的数量。

尽管SIMD指令集已经集成到普通CPU中相当长一段时间了,但是这些扩展功能(SIMD指令集)只能通过在C或者C++代码中通过准汇编语言进行使用。像GNU编译器系列(GCC)这样的现代编译器现在可以将常规的C或C
++源代码转换成向量操作。这个过程被称为自动矢量化(auto-vectorization)。
这篇文章,我们将会介绍必要的软件工具和编程技术,并进一步提供利用GCC自动矢量化功能的示例源代码

在逻辑回归中,以计算z为例,$ z =  w^{T}+b $,你可以用for循环来实现。

硬件


要估算自动矢量化代码可以带来多大的性能提升,首先需要了解您正在使用的CPU的能力以及CPU可以处理多少个Vector。因此,运行以下命令:

$ cat /proc/cpuinfo…
flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat
pse36 clflush dts acpi mmx fxsr sse sse2ss ht tm pbe syscall
nx rdtscp lm constant_tsc arch_perfmon pebs bts rep_good nopl
xtopology nonstop_tsc aperfmperf pni pclmulqdq dtes64 monitor ds_cpl
vmx est tm2ssse3cx16 xtpr pdcmsse4_1 sse4_2x2apic popcnt
xsaveavxlahf_lm ida arat epb xsaveopt pln pts dts tpr_shadow vnmi
flexpriority ept vpid…

在CPU属性中,您将找到有关SIMD功能的信息(以粗体突出显示)。若信息中包含下表列出的字符串,则您的CPU支持(SIMD)功能:

Name Register Size Amount of floats Amount of doubles
mmx (*) 64 bit 0 0
sse2 128 bit 4 2
ssse3 128 bit 4 2
sse4_1 128 bit 4 2
sse4_2 128 bit 4 2
avx 256 bit 8 4
  • 英特尔处理器上的第一个矢量单元只能容纳两个迭代器。

所有Intel和AMD 64位CPU都至少支持sse2指令集。
AVX指令集已于2011年推出,采用Sandy Bridge架构,可在最新一代Mac Book
Pro笔记本电脑上使用。
AVX指令集可以同时对四个双精度浮点值应用数学运算。因此,与未使用矢量单位的版本相比,应用程序可以实现的最大性能增益是4倍。
Vector单位大小的持续增大,反应了是计算机微体系结构演化的新趋势。十年前,计算机的时钟频率几乎每年翻一番。计算机的时钟频率已经达到了极限,制造商转而开始提供越来越大的矢量单位,并且在每个芯片中提供越来越多的核心数量。

但是在python中z可以调用numpy的方法,直接一句$z = np.dot(w,x) +
b$用向量化完成,而且你会发现这个非常快。

哪些代码可以被auto-vectorized ?


为了将计算分布到CPU的矢量单元,编译器必须对源代码的依赖性和相互影响有一个很好的理解。但最重要的是,编译器必须能够检测出那些代码段可以使用SIMD指令进行优化。最简单的情况是对数组执行计算的循环:

for ( unsigned int i = 0; i < ArraySize; ++ i)
{
      c[i] = a[i] * b[i];
}

使用等于或者高于gcc461(在slc5_amd64_gcc461
scram体系结构的服务器中可用)版本的gcc编译器,编译的时候,一定要带上is-ftree-vectorize

注意:从自动矢量化中获益最多的循环是包含数学运算的循环。只是迭代对象集合的循环不会获利,在大多数情况下甚至不能自动矢量化。

可以在这里找到可以自动矢量化的循环结构列表:http://gcc.gnu.org/projects/tree-ssa/vectorization.html
一旦你使用选项-ftree-vectorizer-verbose =
7
编译你的代码,GCC将会给你一个关于你的程序中所有循环的详细报告,以及它们是否已经被自动矢量化了。以下报告是成功对循环进行向量化的结果:

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_get_data_access_cost: inside_cost = 1, outside_cost = 0.

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_get_data_access_cost: inside_cost = 2, outside_cost = 0.

autovect.cpp:66: note: vect_model_store_cost: aligned.

autovect.cpp:66: note: vect_get_data_access_cost: inside_cost = 3, outside_cost = 0.

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_model_load_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: vect_model_load_cost: aligned.

autovect.cpp:66: note: vect_model_load_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: vect_model_simple_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: vect_model_simple_cost: inside_cost = 1, outside_cost = 1 .

autovect.cpp:66: note: vect_model_store_cost: aligned.

autovect.cpp:66: note: vect_model_store_cost: inside_cost = 1, outside_cost = 0 .

autovect.cpp:66: note: Cost model analysis:

Vector inside of loop cost: 5

Vector outside of loop cost: 11

Scalar iteration cost: 5

Scalar outside cost: 0

prologue iterations: 0

epilogue iterations: 2

Calculated minimum iters for profitability: 3

autovect.cpp:66: note:  Profitability threshold = 3

autovect.cpp:66: note: Profitability threshold is 3 loop iterations.

autovect.cpp:66: note: LOOP VECTORIZED.

若一个循环不能被向量化, GCC将会给出原因:

autovect.cpp:133: note: not vectorized: control flow in loop.

上面的输出的含义是:在循环中调用了tostd ::
cout,引入了一个控制流,从而无法对该循环进行向量化。

ng做了个实验,求两个100万长的一维向量的內积,用向量化花了1.5毫秒,而用for循环计算花了400多毫秒。

让GCC自动对你的代码进行向量化


想要让GCC对自己的代码进行自动向量化,需要使用最新的编译器,建议使用至少GCC
4.6.1
版本的编译器,通常来说,版本越新越好。

所以平常记得用向量化,一定要避免使用for循环,你的代码会快很多。

编译器标示


想要打开auto-vectorization,使用标示:

-ftree-vectorize

如果您使用优化级别-O3or进行编译,则隐式启用此选项。
想要得到哪些loop是已经被auto-vectorized和哪些loops没有被成功矢量化以及原因,可以使用选项:

-ftree-vectorizer-verbose=7

关于如何阅读这个输出,请看下面的Real-World代码示例部分。
某些循环结构(例如浮点数的减法)只能在允许编译器改变数学运算顺序的情况下进行矢量化。要做到这一点,需要使用选项:

-ffast-math

如果您使用优化级别-Ofast,则会隐式启用该选项。请注意,fast-math选项会修正修改浮点运算中的错误操作。详情请看http://gcc.gnu.org/onlinedocs/gcc-4.1.1/gcc/Optimize-Options.html
此外,您可以明确指出编译器在生成二进制文件时应使用哪一个SIMD指令集。在编译x86-64架构时,GCC默认使用SSE2指令集。如果要启用AVX指令集,请使用编译器标志:

-mavx

需要注意的是,要运行二进制文件的机器必须支持AVX指令集。您还可以让通过下面的选项,让编译器自己决定使用机器中的哪个指令集:

-march=native

若你想要使用 C++11 的特性,例如:lambda 表达式,一定要确认开启了新的C++
Standard:

-std=gnu++0x

CPU和GPU都有并行化的指令,有时候叫SIMD( single instruction multiple data
)。

启用auto-vectorization最佳实践


如果你使用了这样的内置函数,比如np.function,python的numpy能充分利用并行化去更快的计算。

使用数组结构


想要启用自动向量化,编译器不仅需要能够分析循环迭代,还需要能够分析访问内存的模式。嵌套数组的复杂类型很难或不可能由编译器自动进行矢量化。建议使用简单的c数组orstd
:: arrays(在C ++
11中引入)来保存数据,并允许编译器以连续的方式访问数据。这也将使您的程序能够利用CPU的各种缓存。

如果你需要可变数组,建议使用std :: vectorand
获取指向第一个元素的指针来访问循环中的元素:

std::vectormyData(23);

double * pMyData = &myData.front();

...

pMyData[i] += 5;

数组结构的完整示例可以在下面的“综合代码示例”一节中找到。

 

限制分支


尽力限制for循环内的分支。
GCC能够将一些if语句翻译成向量化的代码,但很有限。这样做时,将评估所有可能的分支,并丢弃未采用分支的值。

2 更多向量化的例子( More Vectorization Examples )

尽量简化代码


如果在for循环中有复杂的计算,您希望GCC对其进行矢量化,可以考虑使用C ++
11中的lambda表达式,将它们分解为更简单的代码。在这里可以找到C
++11新功能的介绍:http://en.wikipedia.org/wiki/Anonymous_function#C.2B.2B

这里是一个lambda函数的例子,用来计算一个double变量的平方。

auto kernel_square =

[] // capture nothing

( double const& a) // take 1 parameter by reference

->double // lambda function returns the square

{

return ( a * a );

};

在循环中调用lambda函数:

for ( unsigned int i = 0; i < size; ++ i)

{

ptr_square[i] = kernel_square( ptr_x[i] ) + kernel_square( ptr_y[i] ) + kernel_square( ptr_z[i] ) ;

}

请注意,此代码将由GCC自动矢量化。对lambda函数的调用不会像常规函数调用那样,造成特别大的开销,因为代码会被完全内联。

另一个更全面的例子是lambda函数中的for循环。这个循环也会被GCC自动矢量化:

// Defining a compute kernel to encapsulate a specific computation

auto kernel_multiply =

[ &cFixedMultiply ] // capture the constant by reference of the scope of the lambda expression

( DataArray const& a, DataArray const& b, DataArray & res ) // take 3 parameters by reference

->void // lambda function returns void

{

// simple loop vectorized

for( unsigned int i = 0; i < a.size(); ++ i)

{

res[i] = a[i] * b[i];

}

};

平时要避免使用for循环,善用python的numpy库中的内置函数。

不要调用外部函数


调用外部函数,如exp(),log()等等,会导致循环不能被向量化。因此,您必须决定循环中的数学运算是否足够,从而将循环分解。这意味着,在第一个循环中,计算调用exp()时,需要使用的参数,并将此结果存储在临时数组中。GCC会自动矢量化这个循环。第二个循环将简单地执行对exp()的调用并存储结果。

如果要调用的函数是由您控制的,请尽量尝试将此函数指定为C ++ 11
lambda表达式。可以在这里找到该特性的介绍:http://en.wikipedia.org/wiki/Anonymous_function#C.2B.2B

比如矩阵A和向量v的內积,可以用np.dot。对一列向量v实施指数运算,可以用np.exp,还有各种np.log,np.abs,np.maxmum(
v, 0)等等。

在循环中使用普通的整数计数器


使用普通的整数计数器来构建for循环,而不是std ::
iterator,因为这些使得GCC难以分析内存访问和循环的属性,如迭代计数。

for (vector::iterator it = y.begin(); it != y.end(); it++)

{

(*it) += 23;

}

变为:

for (unsigned int i = 0; i < y.size(); ++i)

{

y[i] += 23;

}

对于 v**2, 1/v这样的操作也要考虑用np里的函数。

综合代码示例


下面的代码给出了一些关于如何用GCC编写自动矢量化C
++代码的例子。复制并粘贴源代码并使用下面命令进行编译即可:

g++ -Ofast -ftree-vectorizer-verbose=7 -march=native -std=c++11 -o autovect autovect.cpp

请保证使用的gcc版本至少是:GCC 4.7

/*

   Compile with ( use at least gcc 4.7 ):
   g++ -Ofast -ftree-vectorizer-verbose=7 -march=native -std=c++11 -o autovect autovect.cpp

*/

#include <math.h>

#include <string>
#include <iostream>
#include <array>
#include <vector>

// Sturcture-Of-Array to hold coordinates
struct Vector3
{
   std::vector<double> x;
   std::vector<double> y;
   std::vector<double> z;

   // final result of the distance calcualtion
   std::vector<double> distance;

   void add( double _x, double _y, double _z )
   {
      x.push_back( _x );
      y.push_back( _y );
      z.push_back( _z );
      distance.push_back( 0.0f );
   }

};


int main()
{

// Fixed Size Arrays

   typedef std::array<double, 10> DataArray;

   DataArray vect_a = { 0,1,2,3,4,5,6,7,8,9 };
   DataArray vect_b = {0.5,1,2,3,4,5,6,7,8,9 };
   DataArray vect_res_plain = { 0,0,0,0,0,0,0,0,0,0};   
   DataArray vect_res_lambda = { 0,0,0,0,0,0,0,0,0,0};

   constexpr double cFixedMultiply = 23.5f;

   // simple loop vectorized
   // -- auto-vectorized --
   for( unsigned int i = 0; i < vect_a.size(); ++ i)
   {
      vect_res_plain[i] = vect_a[i]  + vect_b[i];
   }

   // Defining a compute kernel to encapsulate a specific computation
   auto kernel_multiply = 
      [ &cFixedMultiply ] // capture the constant by reference of the scope of the lambda expression
      ( DataArray const& a, DataArray const& b, DataArray & res ) // take 3 parameters by reference
      ->void // lambda function returns void
   {
      // simple loop vectorized
      // -- auto-vectorized --
      for( unsigned int i = 0; i < a.size(); ++ i)
      {
         res[i] = a[i] * b[i] * cFixedMultiply;
      }
   };

   // call the lambda function
   // this call is autovectorized
   kernel_multiply ( vect_a, vect_b, vect_res_lambda );


   // This kernel will be called multiple times and performs the quadrature
   auto kernel_square = 
      [] // capture nothing
      ( double const& a) // take 1 parameter by reference
      ->double // lambda function returns the square
   {
      return ( a * a );
   };

   // create struct and fill with dummy values
   Vector3 v3;
   for ( unsigned int i = 0; i < 50 ; ++ i)
   {
      v3.add( i * 1.1, i * 2.2,  i* 3.3 );   
   }

   // store the size in a local variable. This is needed to GCG
   // can estimate the loop iterations
   const unsigned int size = v3.x.size();

   // -- auto-vectorized --
   for ( unsigned int i = 0; i < size; ++ i)
   {
      v3.distance[i] = sqrt( kernel_square( v3.x[i] ) + kernel_square( v3.y[i] ) + kernel_square( v3.z[i] )) ;
   }

   // output the result, so GCC won't optimize the calculations away
   std::cout << std::endl << "Computation on std::array" << std::endl;
   for( unsigned int i = 0; i < vect_a.size(); ++ i)
   {
      std::cout << vect_res_plain[i] << std::endl;
      std::cout << vect_res_lambda[i] << std::endl;
   }

   std::cout << std::endl << "Computation on Structure-of-Array with variable sized std::vectors" << std::endl;
   for( unsigned int i = 0; i < v3.x.size(); ++ i)
   {
      std::cout << "sqrt( " << v3.x[i] << "^2 + " << v3.y[i] << "^2 + " << v3.z[i] << "^2 ) = " 
                << v3.distance[i] << std::endl;
   }

   return 0;
}

图片 4

自动矢量化的实际应用


在这里,您可以找到一部分利用GCC编译器的自动矢量化功能的应用列表:
The VDT mathematical
library

参考:
https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookWritingAutovectorizableCode

图片 5

 

3 向量化逻辑回归( Vectorizing Logistic Regression )

对于逻辑回归的导数计算也应该使用向量化,完全不用for循环。图中给出了向量化的过程。

Z的计算的向量化形式是$z = np.dot(w.T,x) +
b$,其中b在这里是一个实数,python在向量和实数相加时,会自动把实数变成一个相同维度的向量再相加。

 其中w是n * 1的列向量,w.T是1 * n的列向量,X是n * m的矩阵,结果就是1
* m的向量,最后加上1 * m的b向量,得到1 *
m的Z。最后通过sigmoid得到预测值A。

图片 6

 

同时还可以利用向量化计算m个数据的梯度,注意是同时计算。下图左边是for循环的实现,右边是向量化的实现。

这里dz是代价函数对z变量的导数,之前推导过等于预测值减去实际值a – y。

dw是代价函数对w的导数,db是代价函数对b的导数,如果不记得了可以翻看上一节课,逻辑回归的内容。

虽然要尽量使用向量化,但是在进行多次梯度下降的迭代还是要用到for循环,这个不可避免。

图片 7

 

4 python中的广播( python broadcasting)

当你用一个向量加上一个数的时候,python会自动把这个数变成向量再一一相加。

当你用一个m*n的矩阵加(减乘除)上1*n的向量时,python会自动把1*n的向量竖直复制变成m*n再相加。

当你用一个m*n的矩阵加上m*1的向量时,python会自动把m*1的向量水平复制变成m*n再相加。

这是实现神经网络时主要用到的广播,更详细的可以查看numpy文档搜索broadcasting。

对于numpy中的一些用法需要了解,可以帮助你更高效地用矩阵运算来提升程序效率,ng在本节还举了求百分比的例子。

$A.sum(axis=0)$代表竖直求和,如果axis = 1就是水平求和。

图片 8

图片 9

 

5 python / numpy中的向量说明( A note on python/numpy vectors )

numpy和广播使我们可以用一行代码完成很多运算。

但有时可能会引入非常细微的错误,非常奇怪的bug,如果你不熟悉所有的复杂的广播运作方式。

比如你觉得一个行向量和列向量相加应该会报错,但是并不会,而且也不是简单的一一相加。

python这些奇怪的效果有其内在逻辑,如果不熟悉python,你可能会写出奇怪的难以调试的bug。

ng的建议,在实现神经网络的时候不要使用shape为(n,)这样的变量,要用(n,1)。

比如a 的 shape是(5, ) ,当你计算$np.dot(a,
a.T)$的时候得到的是一个实数,a和a的转置,它们的shape都是(5, )。

如果a 的 shape是(5, 1),你计算$np.dot(a,
a.T)$的时候得到的就是一个5*5的矩阵。a的shape是( 5, 1),而a.T的shape是(
1, 5 )。

a.shape = (5,
)这是一个秩为1的数组,不是行向量也不是列向量。很多学生出现难以调试的bug都来自秩为1数组。

另外你在代码中做了很多事情后可能不记得或者不确定a是怎样的时候,用$assert(
a.shape == (5,1) )$来检查你的矩阵的维度。

如果你得到了(5,) 你可以把它reshape成(5, 1)或(1,
5),reshape是很快的O(1)复杂度,所以放心大胆的用它,不用担心。

图片 10