c++ - Cuda:最小二乘求解，速度较差

Question

最近，我使用 Cuda 编写了一个名为“正交匹配追踪”的算法。在我丑陋的 Cuda 代码中，整个迭代需要 60 秒，而 Eigen lib 只需要 3 秒......

在我的代码中，矩阵 A 是 [640,1024]，y 是 [640,1]，在每一步中，我从 A 中选择一些 vector 来组成一个名为 A_temp [640,itera] 的新矩阵，iter=1:500 .我在 cpu 中新建了一个数组 MaxDex_Host[] 来告诉选择哪一列。

我想使用最小二乘法从 A_temp*x_temp=y 获取 x_temp[itera,1]，我使用 cula API“culaDeviceSgels”和 cublas 矩阵 vector 乘法 API。

所以 culaDeviceSgels 会调用 500 次，我认为这会比 Eigen lib 的 QR.Sovler 更快。

我检查了 Nisight 性能分析，我发现 custreamdestory 需要很长时间。我在迭代前初始化 cublas，并在得到结果后销毁它。所以我想知道什么是 custreamdestory ，它与 cublasdestory 有什么不同？

主要问题是 memcpy 和函数 'gemm_kernel1x1val' 。我认为这个函数来自“culaDeviceSgels”

while(itera<500): 我使用 cublasSgemv 和 cublasIsamax 来获取 MaxDex_Host[itera] ，然后

        MaxDex_Host[itera]=pos;
    itera++; 
    float* A_temp_cpu=new float[M*itera]; // matrix all in col-major
    for (int j=0;j<itera;j++) // to  get A_temp [M,itera] , the MaxDex_Host[] shows the positon of which column of A to chose , 
    {
        for (int i=0;i<M;i++) //M=640 , and A is 640*1024 ,itera is add 1 each step
        {
            A_temp_cpu[j*M+i]=A[MaxDex_Host[j]*M+i];
        }
    }
          // I must allocate one more array because culaDeviceSgels will decompose the one input Array ,  and I want to use A_temp after least-square solving.
    float* A_temp_gpu;
    float* A_temp2_gpu;  
    cudaMalloc((void**)&A_temp_gpu,Size_float*M*itera);
    cudaMalloc((void**)&A_temp2_gpu,Size_float*M*itera);
    cudaMemcpy(A_temp_gpu,A_temp_cpu,Size_float*M*itera,cudaMemcpyHostToDevice);
    cudaMemcpy(A_temp2_gpu,A_temp_gpu,Size_float*M*itera,cudaMemcpyDeviceToDevice);
    culaDeviceSgels('N',M,itera,1,A_temp_gpu,M,y_Gpu_temp,M);// the x_temp I want is in y_Gpu_temp's return value ,  stored in the y_Gpu_temp[0]——y_Gpu_temp[itera-1]
     float* x_temp;
    cudaMalloc((void**)&x_temp,Size_float*itera);
    cudaMemcpy(x_temp,y_Gpu_temp,Size_float*itera,cudaMemcpyDeviceToDevice);

Cuda的内存管理似乎太复杂了，有没有其他方便的方法来解决最小二乘问题？

Answer 1

我认为 custreamdestory 和 gemm_kernel1x1val 是由您正在使用的 API 内部调用的，因此与它们没有太大关系。

为了改进您的代码，我建议您执行以下操作。

1. 您可以通过保留矩阵A 的设备拷贝来摆脱A_temp_cpu。然后，您可以通过内核分配将 A 的行复制到 A_temp_gpu 和 A_temp2_gpu 的行中。这将避免执行前两个 cudaMemcpy。
2. 您可以在 while 循环之外预分配 A_temp_gpu 和 A_temp2_gpu，方法是使用 itera 的最大可能值迭代。这将避免循环中的前两个 cudaMalloc。这同样适用于 x_temp。
3. 据我所知，culaDeviceSgels 求解线性方程组。我认为您也可以仅使用 cuBLAS API 来做同样的事情。例如，您可以先通过 cublasDgetrfBatched() 执行 LU 因式分解，然后使用两次 cublasStrsv() 来求解两个产生的线性系统。您可能希望了解此解决方案是否会导致更快的算法。