algorithm - renglón - hallar una base para el espacio columna de una matriz
Cálculo del espacio nulo de una matriz lo más rápido posible. (7)
Necesito calcular el espacio nulo de varios miles de matrices pequeñas (8x9, no 4x3 como escribí anteriormente) en paralelo (CUDA). Todas las referencias apuntan a SVD, pero el algoritmo en las recetas numéricas parece muy caro y me da muchas otras cosas además del espacio nulo que realmente no necesito. ¿La eliminación gaussiana no es realmente una opción? ¿Existen otros métodos de uso común?
"parece muy caro": ¿qué datos tiene que lo respalde?
Tal vez Block Lanczos sea la respuesta que buscas.
O tal vez this .
Tanto JAMA como Apache Commons Math tienen implementaciones SVD en Java. ¿Por qué no tomarlas y probarlas? Obtenga algunos datos reales para su caso en lugar de impresiones. No le costará mucho, ya que el código ya está escrito y probado.
Creo que lo más importante para CUDA es encontrar un algoritmo que no dependa de la bifurcación condicional (que es bastante lenta en el hardware de gráficos). Simple si las declaraciones que pueden optimizarse en la asignación condicional son mucho mejores (o puede usar el operador?:).
Si es necesario, deberías poder hacer algún tipo de giro usando la asignación condicional. En realidad, podría ser más difícil determinar cómo almacenar su resultado: si su matriz es deficiente en rango, ¿qué desea que haga su programa CUDA al respecto?
Si asume que su matriz 4x3 no es realmente deficiente en rango, puede encontrar su vector (simple) de espacio nulo sin ningún condicional: la matriz es lo suficientemente pequeña como para que pueda usar la regla de Cramer de manera eficiente.
En realidad, como no te importa la escala de tu vector nulo, no tienes que dividir por el determinante, solo puedes tomar los determinantes de los menores:
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|
Tenga en cuenta que estos determinantes 3x3 son solo productos triples; puede guardar el cálculo reutilizando los productos cruzados.
En las respuestas anteriores, ya se ha señalado cómo se puede calcular el espacio nulo de una matriz utilizando el enfoque QR o SVD. Se debe preferir SVD cuando se requiere precisión, vea también Espacio nulo de una matriz densa rectangular .
A partir de febrero de 2015, CUDA 7 (ahora en versión candidata) hace disponible la SVD a través de su nueva biblioteca cuSOLVER. A continuación informo un ejemplo sobre cómo usar la SVD de cuSOLVER para calcular el espacio nulo de una matriz.
Tenga en cuenta que el problema en el que se está enfocando se relaciona con el cálculo de varias matrices pequeñas, por lo que debe adaptar el ejemplo que se proporciona a continuación utilizando las secuencias para que tenga sentido para su caso. Para asociar un flujo a cada tarea puede utilizar
cudaStreamCreate()
y
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;
}
Utilidades.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
Utilidades.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://.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__); }
La eliminación gaussiana es bastante rápida para matrices 4x3. IIRC He hecho unos 5 millones por segundo con Java sin paralelismo. Con un problema tan pequeño, lo mejor que puede hacer es codificar la rutina (reducir filas, etc.); de lo contrario, perderá la mayor parte del tiempo poniendo los datos en el formato correcto para la rutina externa.
Me pregunté si las matrices están relacionadas en lugar de ser aleatorias, por lo que los espacios nulos que está buscando pueden considerarse como tangentes unidimensionales a una curva en el espacio N (N = 9). Si es así, puede acelerar las cosas utilizando el método de Newton para resolver instancias sucesivas del sistema de ecuaciones cuadráticas Ax = 0, | x | ^ 2 = 1, comenzando desde un vector de espacio nulo anterior. El método de Newton utiliza los primeros derivados para converger en una solución, y por lo tanto usaría la eliminación gaussiana para resolver sistemas 9x9. El uso de esta técnica requeriría poder realizar pequeños pasos de matriz a matriz, por ejemplo, variando un parámetro.
Entonces, la idea es que inicie utilizando SVD en la primera matriz, pero luego pase de una matriz a otra, utilizando el vector de espacio nulo de uno como punto de partida para la iteración de la siguiente. Necesita una o dos iteraciones para obtener la convergencia. Si no obtiene convegencia, use SVD para reiniciar. Si esta situación es la que tiene, es mucho más rápido que comenzar de nuevo en cada matriz.
Utilicé esto hace mucho tiempo para mapear contornos en las soluciones de conjuntos de 50 x 50 ecuaciones cuadráticas asociadas con el comportamiento de los sistemas de energía eléctrica.
No creo que el método propuesto anteriormente le dé todo el espacio nulo. Para recapitular: "A = QR, donde Q = [Q1 Q2], y Q1 es m-by-n y Q2 es m-by- (mn). Entonces las columnas de Q2 forman el espacio nulo de A ^ T".
De hecho, esto solo puede dar un subespacio del espacio nulo. El contraejemplo simple es cuando A = 0, en cuyo caso el espacio nulo de A ^ T es el valor de R ^ m.
Por lo tanto, es necesario comprobar R también. De acuerdo con mi experiencia con Matlab, si una fila de R es 0, la columna correspondiente en Q también debería ser una base del espacio nulo de A ^ T. Claramente, esta observación es heurística y depende del algoritmo particular utilizado para la descomposición QR.
Para responder a su pregunta directamente ... ¡sí! Descomposición QR!
Sea A una matriz de m por n con rango n. La descomposición QR encuentra la matriz oronal normal m por m Q y la matriz triangular superior m n por r tal que A = QR. Si definimos Q = [Q1 Q2], donde Q1 es m-by-n y Q2 es m-by- (mn), entonces las columnas de Q2 forman el espacio nulo de A ^ T.
La descomposición QR se calcula mediante Gram-Schmidt, rotaciones de Givens o reflexiones de los propietarios. Tienen diferentes propiedades de estabilidad y conteos de operación.
Tienes razón: SVD es caro! No puedo hablar por lo que usa el estado de la técnica, pero cuando escucho "calcular espacio nulo" (EDITAR: de una manera que es simple de entender para mí), creo que QR.