6

我在另一个线程上提出了一个类似的问题，但后来我专注于如何使用 OpenCV。未能达到我最初想要的，我将在这里确切地问我想要什么。

我有两个矩阵。矩阵 a 是 2782x128，矩阵 b 是 4000x128，都是 unsigned char 值。这些值存储在单个数组中。对于 a 中的每个向量，我需要 b 中具有最近欧几里德距离的向量的索引。

好的，现在我的代码来实现这一点：

#include <windows.h>
#include <stdlib.h>
#include <stdio.h>
#include <cstdio>
#include <math.h>
#include <time.h>
#include <sys/timeb.h>
#include <iostream>
#include <fstream>
#include "main.h"

using namespace std;

void main(int argc, char* argv[])
{
    int a_size;
    unsigned char* a = NULL;
    read_matrix(&a, a_size,"matrixa");
    int b_size;
    unsigned char* b = NULL;
    read_matrix(&b, b_size,"matrixb");

    LARGE_INTEGER liStart;
    LARGE_INTEGER liEnd;
    LARGE_INTEGER liPerfFreq;
    QueryPerformanceFrequency( &liPerfFreq );
    QueryPerformanceCounter( &liStart );

    int* indexes = NULL;
    min_distance_loop(&indexes, b, b_size, a, a_size);

    QueryPerformanceCounter( &liEnd );

    cout << "loop time: " << (liEnd.QuadPart - liStart.QuadPart) / long double(liPerfFreq.QuadPart) << "s." << endl;

    if (a)
    delete[]a;
if (b)
    delete[]b;
if (indexes)
    delete[]indexes;
    return;
}

void read_matrix(unsigned char** matrix, int& matrix_size, char* matrixPath)
{
    ofstream myfile;
    float f;
    FILE * pFile;
    pFile = fopen (matrixPath,"r");
    fscanf (pFile, "%d", &matrix_size);
    *matrix = new unsigned char[matrix_size*128];

    for (int i=0; i<matrix_size*128; ++i)
    {
        unsigned int matPtr;
        fscanf (pFile, "%u", &matPtr);
        matrix[i]=(unsigned char)matPtr;
    }
    fclose (pFile);
}

void min_distance_loop(int** indexes, unsigned char* b, int b_size, unsigned char* a, int a_size)
{
    const int descrSize = 128;

    *indexes = (int*)malloc(a_size*sizeof(int));
    int dataIndex=0;
    int vocIndex=0;
    int min_distance;
    int distance;
    int multiply;

    unsigned char* dataPtr;
    unsigned char* vocPtr;
    for (int i=0; i<a_size; ++i)
    {
        min_distance = LONG_MAX;
        for (int j=0; j<b_size; ++j)
        {
            distance=0;
            dataPtr = &a[dataIndex];
            vocPtr = &b[vocIndex];

            for (int k=0; k<descrSize; ++k)
            {
                multiply = *dataPtr++-*vocPtr++;
                distance += multiply*multiply;
                // If the distance is greater than the previously calculated, exit
                if (distance>min_distance)
                    break;
            }

            // if distance smaller
            if (distance<min_distance)
            {
                min_distance = distance;
                (*indexes)[i] = j;
            }
            vocIndex+=descrSize;
        }
        dataIndex+=descrSize;
        vocIndex=0;
    }
}

并附有样本矩阵的文件。

矩阵 a 矩阵b

我使用windows.h只是为了计算消耗时间，所以如果你想在windows以外的其他平台测试代码，只需更改windows.h头文件并改变计算消耗时间的方式。

我电脑中的这段代码大约是 0.5 秒。问题是我在 Matlab 中有另一个代码可以在 0.05 秒内完成同样的事情。在我的实验中，我每秒收到几个矩阵，比如矩阵 a，所以 0.5 秒太多了。

现在计算这个的matlab代码：

aa=sum(a.*a,2); bb=sum(b.*b,2); ab=a*b'; 
d = sqrt(abs(repmat(aa,[1 size(bb,1)]) + repmat(bb',[size(aa,1) 1]) - 2*ab));
[minz index]=min(d,[],2);

好的。Matlab 代码使用的是 (xa)^2 = x^2 + a^2 - 2ab。

所以我的下一个尝试是做同样的事情。我删除了自己的代码以进行相同的计算，但大约是 1.2 秒。

然后，我尝试使用不同的外部库。第一次尝试是 Eigen：

const int descrSize = 128;
MatrixXi a(a_size, descrSize);
MatrixXi b(b_size, descrSize);
MatrixXi ab(a_size, b_size);

unsigned char* dataPtr = matrixa;
for (int i=0; i<nframes; ++i)
{
    for (int j=0; j<descrSize; ++j)
    {
        a(i,j)=(int)*dataPtr++;
    }
}
unsigned char* vocPtr = matrixb;
for (int i=0; i<vocabulary_size; ++i)
{
    for (int j=0; j<descrSize; ++j)
    {
        b(i,j)=(int)*vocPtr ++;
    }
}
ab = a*b.transpose();
a.cwiseProduct(a);
b.cwiseProduct(b);
MatrixXi aa = a.rowwise().sum();
MatrixXi bb = b.rowwise().sum();
MatrixXi d = (aa.replicate(1,vocabulary_size) + bb.transpose().replicate(nframes,1) - 2*ab).cwiseAbs2();

int* index = NULL;
index = (int*)malloc(nframes*sizeof(int));
for (int i=0; i<nframes; ++i)
{
    d.row(i).minCoeff(&index[i]);
}

这个特征码的成本大约为 1.2，仅用于以下行： ab = a*b.transpose();

还使用了使用 opencv 的类似代码，并且 ab = a*b.transpose(); 的成本 为 0.65 秒。

所以，matlab 能够这么快地做同样的事情，而我在 C++ 中却做不到，这真的很烦人！当然，能够运行我的实验会很棒，但我认为缺乏知识才是真正让我烦恼的地方。我怎样才能达到至少与 Matlab 相同的性能？欢迎任何形式的解决方案。我的意思是，任何外部库（如果可能的话免费）、循环展开的东西、模板的东西、SSE 指令（我知道它们存在）、缓存的东西。正如我所说，我的主要目的是增加我的知识，以便能够以更快的性能编写这样的想法。

提前致谢

编辑：David Hammen 建议的更多代码。在进行任何计算之前，我将数组转换为 int 。这是代码：

void min_distance_loop(int** indexes, unsigned char* b, int b_size, unsigned char* a, int a_size)
{
    const int descrSize = 128;

    int* a_int;
    int* b_int;

    LARGE_INTEGER liStart;
    LARGE_INTEGER liEnd;
    LARGE_INTEGER liPerfFreq;
    QueryPerformanceFrequency( &liPerfFreq );
    QueryPerformanceCounter( &liStart );

    a_int = (int*)malloc(a_size*descrSize*sizeof(int));
    b_int = (int*)malloc(b_size*descrSize*sizeof(int));

    for(int i=0; i<descrSize*a_size; ++i)
        a_int[i]=(int)a[i];
    for(int i=0; i<descrSize*b_size; ++i)
        b_int[i]=(int)b[i];

    QueryPerformanceCounter( &liEnd );

    cout << "Casting time: " << (liEnd.QuadPart - liStart.QuadPart) / long double(liPerfFreq.QuadPart) << "s." << endl;

    *indexes = (int*)malloc(a_size*sizeof(int));
    int dataIndex=0;
    int vocIndex=0;
    int min_distance;
    int distance;
    int multiply;

    /*unsigned char* dataPtr;
    unsigned char* vocPtr;*/
    int* dataPtr;
    int* vocPtr;
    for (int i=0; i<a_size; ++i)
    {
        min_distance = LONG_MAX;
        for (int j=0; j<b_size; ++j)
        {
            distance=0;
            dataPtr = &a_int[dataIndex];
            vocPtr = &b_int[vocIndex];

            for (int k=0; k<descrSize; ++k)
            {
                multiply = *dataPtr++-*vocPtr++;
                distance += multiply*multiply;
                // If the distance is greater than the previously calculated, exit
                if (distance>min_distance)
                    break;
            }

            // if distance smaller
            if (distance<min_distance)
            {
                min_distance = distance;
                (*indexes)[i] = j;
            }
            vocIndex+=descrSize;
        }
        dataIndex+=descrSize;
        vocIndex=0;
    }
}

整个过程现在是 0.6，开始的铸造循环是 0.001 秒。也许我做错了什么？

EDIT2：关于本征的任何信息？当我寻找外部库时，他们总是谈论 Eigen 及其速度。我做错了什么？这是一个使用 Eigen 的简单代码，表明它不是那么快。也许我缺少一些配置或一些标志，或者......

MatrixXd A = MatrixXd::Random(1000, 1000);
MatrixXd B = MatrixXd::Random(1000, 500);
MatrixXd X;

这段代码大约是 0.9 秒。

4

3 回答 3

3

正如您所观察到的，您的代码由代表大约 2.8e9 算术运算的矩阵乘积支配。Yopu 说 Matlab（或者更确切地说是高度优化的 MKL）在大约 0.05 秒内计算它。这代表 57 GFLOPS 的速率，表明它不仅使用矢量化，还使用多线程。-fopenmp使用 Eigen，您可以通过启用 OpenMP（使用 gcc）进行编译来启用多线程。在我 5 年前的计算机（2.66Ghz Core2）上，使用浮点数和 4 个线程，您的产品大约需要 0.053 秒，而没有 OpenMP 则需要 0.16 秒，所以您的编译标志肯定有问题。总而言之，要充分利用 Eigen：

  • 64位模式编译
  • 使用浮点数（由于矢量化，双打速度慢了一倍）
  • 启用 OpenMP
  • 如果你的 CPU 有超线程，那么要么禁用它，要么将OMP_NUM_THREADS环境变量定义为物理内核的数量（这很重要，否则性能会很差！）
  • 如果您有其他任务正在运行，最好减少OMP_NUM_THREADS到nb_cores-1
  • 尽可能使用最新的编译器，GCC、clang 和 ICC 最好，MSVC 通常较慢。
于 2013-08-22T09:50:29.467 回答
2

在你的 C++ 代码中肯定会伤害你的一件事是它有大量的 char 到 int 的转换。装船时，我的意思是最多 2*2782*4000*128 char 到 int 的转换。那些char转换int很慢，非常慢。

您可以通过分配一对int数组（一个 2782*128 和另一个 4000*128）将其减少到 (2782+4000)*128 这样的转换，以包含您的char* a和char* b数组的转换为整数的内容。使用这些int*数组而不是您的char*数组。

另一个问题可能是您对intvs的使用long。我不在 Windows 上工作，所以这可能不适用。在我工作的机器上，int是 32 位，long现在是 64 位。32 位绰绰有余，因为 255*255*128 < 256*256*128 = 2 23。

这显然不是问题。

引人注目的是，有问题的代码并未计算 Matlab 代码正在创建的那个巨大的 2728 x 4000 数组。更引人注目的是，Matlab 最有可能使用双精度而不是整数来执行此操作——而且它仍然在 C/C++ 代码中脱颖而出。

一个大问题是缓存。那个 4000*128 数组对于 1 级缓存来说太大了，你正在迭代那个大数组 2782 次。您的代码在内存上等待太多了。要克服这个问题，请使用较小的b数组块，以便您的代码尽可能长时间地使用 1 级缓存。

另一个问题是优化if (distance>min_distance) break;。我怀疑这实际上是一种优化。在if最里面的循环中进行测试通常是个坏主意。尽可能快地冲破那个内积。除了浪费计算之外，摆脱这个测试并没有什么坏处。有时，如果这样做可以删除最内层循环中的分支，则最好进行明显不需要的计算。这是其中一种情况。您可能只需消除此测试即可解决您的问题。尝试这样做。

回到缓存问题，您需要摆脱这个分支，以便可以将a和b矩阵上的操作拆分成更小的块，一次不超过 256 行的块。这就是两个现代 Intel 芯片的 L1 缓存之一中有多少行 128 个无符号字符。由于 250 除以 4000，因此请考虑在逻辑上将该b矩阵拆分为 16 个块。您可能希望形成 2872 x 4000 的大内积数组，但要分小块进行。您可以将其if (distance>min_distance) break;重新添加，但在块级别而不是逐字节级别添加。

您应该能够击败 Matlab，因为它几乎可以肯定是使用双精度数，但您可以使用无符号字符和整数。

于 2012-09-26T09:22:50.693 回答
1

矩阵乘法通常对两个矩阵之一使用最差的缓存访问模式，解决方案是转置其中一个矩阵并使用专门的乘法算法来处理以这种方式存储的数据。

您的矩阵已转置存储。通过将其转换为正常顺序，然后使用正常矩阵相乘，您绝对会扼杀性能。

编写您自己的矩阵乘法循环，将索引的顺序反转到第二个矩阵（具有转置它的效果，而实际上没有移动任何东西并破坏缓存行为）。并通过您的编译器启用自动矢量化的任何选项。

于 2012-09-26T14:18:53.377 回答