Вычисление нулевого пространства матрицы как можно быстрее

Ответ 1

Чтобы ответить на ваш вопрос напрямую... да! QR-разложение!

Пусть A - матрица m-на-n с ранга n. QR-декомпозиция находит ортонормированную m-на-матрицу Q и верхнюю треугольную m-на-матрицу R такую, что A = QR. Если мы определим Q = [Q1 Q2], где Q1 - m-by-n, а Q2 - m-by-(m-n), то столбцы Q2 образуют пустое пространство A ^ T.

QR-декомпозиция вычисляется либо с помощью грамм-Шмидта, Гивенса, либо отражений Домахолдера. Они имеют разные свойства стабильности и операции.

Вы правы: SVD дорого! Я не могу говорить о том, что использует новейшие вещи, но когда я слышу "вычислить нулевое пространство" (EDIT: таким образом, что это просто для меня, чтобы понять), я думаю, что QR.

Ответ 2

Гауссово исключение достаточно быстро для матриц 4x3. IIRC Я сделал около 5 миллионов в секунду с Java без parallelism. С такой небольшой проблемой лучше всего скопировать процедуру (сокращение строк и т.д.); иначе вы будете тратить большую часть времени на то, чтобы данные были в правильном формате для внешней процедуры.

Ответ 3

Я не думаю, что предложенный выше метод всегда дает полное пустое пространство. Напомним: "A = QR, где Q = [Q1 Q2], а Q1 - m-by-n, а Q2 - m-by-(mn). Тогда столбцы Q2 образуют пустое пространство A ^ T."

Действительно, это может дать только подпространство нулевого пространства. Простой встречный пример - это когда A = 0, и в этом случае нулевое пространство A ^ T является целым R ^ m.

Следовательно, необходимо также проверить R. Основываясь на моем опыте с Matlab, если строка R прямая 0, то соответствующий столбец в Q также должен быть базой нулевого пространства A ^ T. Очевидно, что это наблюдение эвристично и зависит от конкретного алгоритма, используемого для QR-декомпозиции.

Ответ 4

"кажется очень дорогим" - какие данные у вас есть, которые поддерживают это?

Возможно, Block Lanczos - ответ, который вы ищете.

Или, может быть, this.

Оба JAMA и Apache Commons Math имеют реализации SVD в Java. Почему бы не взять их и не попробовать? Получите реальные данные для вашего случая, а не впечатления. Это не будет стоить вам многого, так как код уже написан и протестирован.

Ответ 5

Я думаю, что самое главное для 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 детерминанты являются всего лишь тройными продуктами; вы можете сохранить вычисления, повторно используя перекрестные продукты.

Ответ 6

Я задавался вопросом, связаны ли матрицы, а не просто случайны, так что нулевые пробелы, которые вы ищете, можно рассматривать как одномерные касательные к кривой в N-пространстве (N = 9). Если это так, вы можете ускорить процесс, используя метод Ньютона для решения последовательных экземпляров системы квадратичных уравнений Ax = 0, | x | ^ 2 = 1, начиная с предыдущего нулевого пространства-вектора. Метод Ньютона использует первые производные, чтобы сходиться к решению, и поэтому использовал бы гауссово исключение для решения систем 9х9. Использование этого метода потребовало бы, чтобы вы могли делать небольшие шаги от матрицы к матрице, например, изменяя параметр.

Итак, идея состоит в том, что вы инициализируете использование SVD в первой матрице, но после этого вы переходите от матрицы к матрице, используя нулевой вектор пространства одного в качестве отправной точки для итерации для следующей. Для достижения конвергенции вам понадобится одна или две итерации. Если вы не получаете выпуклость, вы используете SVD для перезапуска. Если эта ситуация такова, что у вас есть, она намного быстрее, чем начинать свежую на каждой матрице.

Я использовал это давным-давно, чтобы отобразить контуры в решениях множеств 50 x 50 квадратичных уравнений, связанных с поведением электрических систем.

Ответ 7

В вышеприведенных ответах уже указано, как нулевое пространство матрицы можно вычислить с использованием подхода QR или SVD. SVD следует выбирать, когда требуется точность, см. Также Нулевое пространство прямоугольной плотной матрицы.

По состоянию на февраль 2015 года CUDA 7 (теперь в кандидате на выпуск) делает SVD доступным через свою новую библиотеку cuSOLVER. Ниже я расскажу о том, как использовать 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/info/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__); }