我需要并行计算数千个小矩阵（8x9，而不是我之前写的 4x3）的零空间（CUDA）。所有引用都指向 SVD，但数值配方中的算法似乎非常昂贵，并且除了我并不真正需要的空空间之外，还给了我很多东西。高斯消除真的不是一种选择吗？还有其他常用的方法吗？
【问题讨论】：
标签： algorithm math matrix cuda
“看起来很贵”——你有什么数据支持这一点？
也许Block Lanczos 是您寻求的答案。
或者this。
JAMA 和 Apache Commons Math 都在 Java 中实现了 SVD。为什么不带上这些并尝试一下呢？为您的案例获取一些真实数据，而不是展示次数。它不会花费您太多，因为代码已经编写和测试。
【讨论】：
高斯消除对于 4x3 矩阵来说非常快。 IIRC 在没有并行性的情况下，我每秒使用 Java 完成了大约 500 万次。对于这么小的问题，最好的办法是自己编写例程（行缩减等）；否则您将浪费大部分时间将数据转换为外部例程的正确格式。
【讨论】：
我认为 CUDA 最重要的是找到一种不依赖于条件分支的算法（这在图形硬件上相当慢）。可以优化为条件赋值的简单 if 语句要好得多（或者您可以使用 ?: 运算符）。
如有必要，您应该能够使用条件赋值进行某种形式的旋转。实际上可能更难确定如何存储您的结果：如果您的矩阵秩不足，您希望您的 CUDA 程序如何处理它？
如果您假设您的 4x3 矩阵实际上不是秩不足的，您可以在没有任何条件的情况下找到您的（单个）零空间向量：矩阵足够小，您可以有效地使用 Cramer 规则。
实际上，由于您实际上并不关心空向量的规模，因此您不必除以行列式 - 您可以只取未成年人的行列式：
x1 x2 x3
M = y1 y2 y3
z1 z2 z3
w1 w2 w3
|y1 y2 y3| |x1 x2 x3| |x1 x2 x3| |x1 x2 x3|
-> x0 = |z1 z2 z3| y0 = -|z1 z2 z3| z0 = |y1 y2 y3| w0 = -|y1 y2 y3|
|w1 w2 w3| |w1 w2 w3| |w1 w2 w3| |z1 z2 z3|
请注意，这些 3x3 行列式只是三倍积；您可以通过重用叉积来节省计算量。
【讨论】：
我想知道矩阵是否相关，而不仅仅是随机的，因此您正在寻找的零空间可以被认为是 N 空间 (N = 9) 中曲线的一维切线。如果是这样，您可以通过使用牛顿法从先前的零空间向量开始求解二次方程组 Ax = 0, |x|^2 = 1 的连续实例来加快速度。牛顿法使用一阶导数收敛到一个解，因此将使用高斯消元法求解 9x9 系统。使用这种技术需要你能够通过改变一个参数来从一个矩阵到另一个矩阵进行小步骤。
所以想法是您在第一个矩阵上使用 SVD 进行初始化，但此后您从一个矩阵到另一个矩阵，使用一个的零空间向量作为下一个迭代的起点。您需要一到两次迭代才能收敛。如果您没有获得收敛，则使用 SVD 重新启动。如果您遇到这种情况，那么它比在每个矩阵上重新开始要快得多。
我很久以前就用它来绘制与电力系统行为相关的 50 x 50 二次方程组解中的等高线。
【讨论】：
直接回答您的问题...是的！二维码分解！
令 A 是一个秩为 n 的 m×n 矩阵。 QR 分解找到正交 m×m 矩阵 Q 和上三角 m×n 矩阵 R 使得 A = QR。如果我们定义 Q = [Q1 Q2]，其中 Q1 是 m×n，Q2 是 m×(m-n)，那么 Q2 的列就形成了 A^T 的零空间。
QR 分解通过 Gram-Schmidt、Givens 旋转或 Householder 反射计算。它们具有不同的稳定性属性和操作次数。
你是对的：SVD 很昂贵！我不能说最先进的东西使用什么，但是当我听到“计算零空间”（编辑：以一种我很容易理解的方式）时，我认为是 QR。
【讨论】：
我不认为上面提出的方法总是给出整个空空间。回顾一下：“A = QR，其中 Q = [Q1 Q2]，Q1 是 m×n，Q2 是 m×(mn)。然后 Q2 的列形成 A^T 的零空间。”
确实，这可能只给出零空间的一个子空间。简单的反例是当 A=0 时，在这种情况下 A^T 的零空间是整个 R^m。
因此，也需要检查 R。根据我对Matlab的经验，如果R的一行是直0，那么Q中的对应列也应该是A^T的零空间的基础。显然，这种观察是启发式的，并且取决于用于 QR 分解的特定算法。
【讨论】：
在上面的答案中，已经指出如何使用 QR 或 SVD 方法计算矩阵的零空间。当需要准确性时，应首选 SVD，另请参阅 Null-space of a rectangular dense matrix。
截至 2015 年 2 月，CUDA 7（现在处于候选版本中）通过其新的 cuSOLVER 库提供 SVD。下面我举一个例子，说明如何使用 cuSOLVER 的 SVD 计算矩阵的零空间。
请注意，您关注的问题涉及几个小矩阵的计算，因此您应该通过使用流来调整我在下面提供的示例，以便对您的情况有意义。要将流与您可以使用的每个任务相关联
cudaStreamCreate()
和
cusolverDnSetStream()
kernel.cu
#include "cuda_runtime.h"
#include "device_launch_paraMeters.h"
#include<iostream>
#include<iomanip>
#include<stdlib.h>
#include<stdio.h>
#include<assert.h>
#include<math.h>
#include <cusolverDn.h>
#include <cuda_runtime_api.h>
#include "Utilities.cuh"
/********/
/* MAIN */
/********/
int main(){
// --- gesvd only supports Nrows >= Ncols
// --- column major memory ordering
const int Nrows = 7;
const int Ncols = 5;
// --- cuSOLVE input/output parameters/arrays
int work_size = 0;
int *devInfo; gpuErrchk(cudaMalloc(&devInfo, sizeof(int)));
// --- CUDA solver initialization
cusolverDnHandle_t solver_handle;
cusolverDnCreate(&solver_handle);
// --- Singular values threshold
double threshold = 1e-12;
// --- Setting the host, Nrows x Ncols matrix
double *h_A = (double *)malloc(Nrows * Ncols * sizeof(double));
for(int j = 0; j < Nrows; j++)
for(int i = 0; i < Ncols; i++)
h_A[j + i*Nrows] = (i + j*j) * sqrt((double)(i + j));
// --- Setting the device matrix and moving the host matrix to the device
double *d_A; gpuErrchk(cudaMalloc(&d_A, Nrows * Ncols * sizeof(double)));
gpuErrchk(cudaMemcpy(d_A, h_A, Nrows * Ncols * sizeof(double), cudaMemcpyHostToDevice));
// --- host side SVD results space
double *h_U = (double *)malloc(Nrows * Nrows * sizeof(double));
double *h_V = (double *)malloc(Ncols * Ncols * sizeof(double));
double *h_S = (double *)malloc(min(Nrows, Ncols) * sizeof(double));
// --- device side SVD workspace and matrices
double *d_U; gpuErrchk(cudaMalloc(&d_U, Nrows * Nrows * sizeof(double)));
double *d_V; gpuErrchk(cudaMalloc(&d_V, Ncols * Ncols * sizeof(double)));
double *d_S; gpuErrchk(cudaMalloc(&d_S, min(Nrows, Ncols) * sizeof(double)));
// --- CUDA SVD initialization
cusolveSafeCall(cusolverDnDgesvd_bufferSize(solver_handle, Nrows, Ncols, &work_size));
double *work; gpuErrchk(cudaMalloc(&work, work_size * sizeof(double)));
// --- CUDA SVD execution
cusolveSafeCall(cusolverDnDgesvd(solver_handle, 'A', 'A', Nrows, Ncols, d_A, Nrows, d_S, d_U, Nrows, d_V, Ncols, work, work_size, NULL, devInfo));
int devInfo_h = 0; gpuErrchk(cudaMemcpy(&devInfo_h, devInfo, sizeof(int), cudaMemcpyDeviceToHost));
if (devInfo_h != 0) std::cout << "Unsuccessful SVD execution\n\n";
// --- Moving the results from device to host
gpuErrchk(cudaMemcpy(h_S, d_S, min(Nrows, Ncols) * sizeof(double), cudaMemcpyDeviceToHost));
gpuErrchk(cudaMemcpy(h_U, d_U, Nrows * Nrows * sizeof(double), cudaMemcpyDeviceToHost));
gpuErrchk(cudaMemcpy(h_V, d_V, Ncols * Ncols * sizeof(double), cudaMemcpyDeviceToHost));
for(int i = 0; i < min(Nrows, Ncols); i++)
std::cout << "d_S["<<i<<"] = " << std::setprecision(15) << h_S[i] << std::endl;
printf("\n\n");
int count = 0;
bool flag = 0;
while (!flag) {
if (h_S[count] < threshold) flag = 1;
if (count == min(Nrows, Ncols)) flag = 1;
count++;
}
count--;
printf("The null space of A has dimension %i\n\n", min(Ncols, Nrows) - count);
for(int j = count; j < Ncols; j++) {
printf("Basis vector nr. %i\n", j - count);
for(int i = 0; i < Ncols; i++)
std::cout << "d_V["<<i<<"] = " << std::setprecision(15) << h_U[j*Ncols + i] << std::endl;
printf("\n");
}
cusolverDnDestroy(solver_handle);
return 0;
}
Utilities.cuh
#ifndef UTILITIES_CUH
#define UTILITIES_CUH
extern "C" int iDivUp(int, int);
extern "C" void gpuErrchk(cudaError_t);
extern "C" void cusolveSafeCall(cusolverStatus_t);
#endif
Utilities.cu
#include <stdio.h>
#include <assert.h>
#include "cuda_runtime.h"
#include <cuda.h>
#include <cusolverDn.h>
/*******************/
/* iDivUp FUNCTION */
/*******************/
extern "C" int iDivUp(int a, int b){ return ((a % b) != 0) ? (a / b + 1) : (a / b); }
/********************/
/* CUDA ERROR CHECK */
/********************/
// --- Credit to http://stackoverflow.com/questions/14038589/what-is-the-canonical-way-to-check-for-errors-using-the-cuda-runtime-api
void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
if (code != cudaSuccess)
{
fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
if (abort) { exit(code); }
}
}
extern "C" void gpuErrchk(cudaError_t ans) { gpuAssert((ans), __FILE__, __LINE__); }
/**************************/
/* CUSOLVE ERROR CHECKING */
/**************************/
static const char *_cudaGetErrorEnum(cusolverStatus_t error)
{
switch (error)
{
case CUSOLVER_STATUS_SUCCESS:
return "CUSOLVER_SUCCESS";
case CUSOLVER_STATUS_NOT_INITIALIZED:
return "CUSOLVER_STATUS_NOT_INITIALIZED";
case CUSOLVER_STATUS_ALLOC_FAILED:
return "CUSOLVER_STATUS_ALLOC_FAILED";
case CUSOLVER_STATUS_INVALID_VALUE:
return "CUSOLVER_STATUS_INVALID_VALUE";
case CUSOLVER_STATUS_ARCH_MISMATCH:
return "CUSOLVER_STATUS_ARCH_MISMATCH";
case CUSOLVER_STATUS_EXECUTION_FAILED:
return "CUSOLVER_STATUS_EXECUTION_FAILED";
case CUSOLVER_STATUS_INTERNAL_ERROR:
return "CUSOLVER_STATUS_INTERNAL_ERROR";
case CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED:
return "CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED";
}
return "<unknown>";
}
inline void __cusolveSafeCall(cusolverStatus_t err, const char *file, const int line)
{
if(CUSOLVER_STATUS_SUCCESS != err) {
fprintf(stderr, "CUSOLVE error in file '%s', line %d\n %s\nerror %d: %s\nterminating!\n",__FILE__, __LINE__,err, \
_cudaGetErrorEnum(err)); \
cudaDeviceReset(); assert(0); \
}
}
extern "C" void cusolveSafeCall(cusolverStatus_t err) { __cusolveSafeCall(err, __FILE__, __LINE__); }
【讨论】：