Calcul de l'espace nul d'une matrice aussi vite que possible

L'élimination gaussienne est assez rapide pour les matrices 4x3. IIRC j'ai fait environ 5 millions par seconde avec Java sans parallélisme. Avec un tel petit problème, votre meilleur pari est de coder la routine(réduction de ligne etc.) vous-même; sinon vous perdrez la plupart du temps à mettre les données dans le bon format pour la routine externe.

Je ne pense pas que la méthode proposée ci-dessus donne toujours tout l'espace nul. Pour récapituler: "A = QR, OÙ Q = [Q1 Q2], et Q1 Est m-by-n et Q2 Est m-by-(m-n). Alors les colonnes de Q2 forment l'espace nul de A^T."

en Effet, cela peut seulement donner un sous-espace de l'espace nul. Un simple contre-exemple est quand A=0, auquel cas l'espace nul de A^T est l'ensemble de R^M.

par conséquent, il est nécessaire de vérifier R trop. Basé sur mon expérience avec Matlab, si une rangée de R est tout droit 0, puis la colonne correspondante en Q devrait également être une base de l'espace nul de A^T. Cette observation est clairement heuristique et dépend de l'algorithme particulier utilisé pour la décomposition QR.

je pense que la chose la plus importante pour CUDA est de trouver un algorithme qui ne dépend pas de branchement conditionnel (ce qui est assez lent sur le matériel graphique). Simple si les énoncés qui peuvent être optimisés en affectation conditionnelle sont beaucoup mieux (ou vous pouvez utiliser le ?: opérateur.)

si nécessaire, vous devriez être en mesure de faire une certaine forme de pivotement en utilisant assignation conditionnelle. Il pourrait en fait être plus difficile de déterminer comment enregistrer votre résultat: si votre matrice est de rang déficientes, que voulez-vous que votre programme CUDA fasse à ce sujet?

si vous supposez que votre matrice 4x3 n'est pas vraiment déficiente en rang, vous pouvez trouver votre vecteur (simple) null-space sans aucune conditionals du tout: la matrice est assez petite que vous pouvez utiliser la règle de Cramer efficacement.

en fait, puisque vous ne vous souciez pas réellement de l'échelle de votre vecteur nul, vous n'avez pas à diviser par le déterminant -- vous pouvez simplement prendre les déterminants de mineurs:

    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|

notez que ces déterminants 3x3 ne sont que des produits triples; vous pouvez sauver le calcul en réutilisant les produits croisés.

"semble très cher" - quelles données Avez-vous qui soutiennent cela?

Peut-être Bloc De Lanczos est la réponse que vous cherchez.

Ou peut-être .

JAMA et Apache Commons Math ont des implémentations SVD en Java. Pourquoi ne pas prendre ceux et de les essayer? Obtenir certaines données réelles pour votre cas au lieu d'impressions. Cela ne vous coûtera pas beaucoup, puisque le code est déjà écrit et testé.

je me suis demandé si les matrices sont liées plutôt que d'être simplement aléatoires, de sorte que les espaces nuls que vous recherchez peuvent être considérés comme des tangentes de Dimension 1 à une courbe dans L'espace N (N = 9). Si oui, vous pouvez être en mesure d'accélérer les choses en utilisant la méthode de Newton pour résoudre les instances successives du système d'équations quadratiques Ax = 0, |x|^2 = 1, à partir d'un vecteur espace nul précédent. La méthode de Newton utilise d'abord dérivés pour converger vers une solution, et donc utiliserait Gaussian élimination pour résoudre les systèmes 9x9. L'utilisation de cette technique exigerait que vous soyez en mesure de faire de petits pas de matrice à matrice en modifiant par exemple un paramètre.

donc l'idée est que vous initialisez en utilisant SVD sur la première matrice, mais ensuite vous passez de matrice en matrice, en utilisant le vecteur d'espace nul de l'un comme point de départ pour l'itération pour la suivante. Vous avez besoin d'une ou deux itérations pour atteindre la convergence. Si vous n'avez pas de convegence, utilisez SVD pour redémarrer. Si ce situation est ce que vous avez, il est beaucoup plus rapide que de commencer à nouveau sur chaque matrice.

j'ai utilisé ceci il y a longtemps pour cartographier les contours dans les solutions des ensembles de 50 x 50 équations quadratiques associées au comportement des systèmes d'énergie électrique.

dans les anwers ci-dessus, il a déjà été souligné comment l'espace nul d'une matrice peut être calculé en utilisant L'approche QR ou SVD. SVD doit être préféré lorsque la précision est requise, voir aussi Null-l'espace d'un rectangle de matrice dense.

depuis février 2015, CUDA 7 (maintenant en version candidate) rend SVD disponible via sa nouvelle bibliothèque cuSOLVER. Ci-dessous, je donne un exemple sur la façon d'utiliser le SVD de cuSOLVER pour calculer l'espace nul d'un matrice.

soyez conscient que le problème sur lequel vous vous concentrez concerne le calcul de plusieurs petites matrices, donc vous devriez adapter l'exemple que je vous donne ci-dessous en utilisant streams pour donner du sens à votre cas. Pour associer un flux pour chaque tâche, vous pouvez utiliser

cudaStreamCreate()

et

cusolverDnSetStream()

kernel.U

#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;

}

Utilitaires.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

Services Publics.U

#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__); }