Loading...
Least squares finds x minimizing ||Ax - b||₂ for overdetermined systems (more equations than unknowns). Two main approaches: (1) Normal equations: solve A^TAx = A^Tb, (2) QR decomposition: factor A = QR, solve Rx = Q^Tb. QR is more numerically stable but slower. Normal equations are faster but square the condition number.
Use optimized least squares routine with automatic QR factorization.
A^TA is SPD, use Cholesky for 2x speedup when κ is small.
Solve many small systems in parallel.
Normal equations are fast but numerically unstable for ill-conditioned A.
void least_squares_normal(cublasHandle_t blas, float* d_A, float* d_b,
float* d_x, int m, int n) {
float *d_ATA, *d_ATb;
cudaMalloc(&d_ATA, n * n * sizeof(float));
cudaMalloc(&d_ATb, n * sizeof(float));
// A^T A (squares condition number!)
cublasSgemm(blas, CUBLAS_OP_T, CUBLAS_OP_N, n, n, m,
&one, d_A, m, d_A, m, &zero, d_ATA, n);
// A^T b
cublasSgemv(blas, CUBLAS_OP_T, m, n, &one, d_A, m, d_b, 1, &zero, d_ATb, 1);
// Solve A^TA x = A^Tb via Cholesky
cusolverDnSpotrf(solver, CUBLAS_FILL_MODE_LOWER, n, d_ATA, n, ...);
cusolverDnSpotrs(solver, CUBLAS_FILL_MODE_LOWER, n, 1, d_ATA, n, d_ATb, n, ...);
cudaMemcpy(d_x, d_ATb, n * sizeof(float), D2D);
}QR-based least squares is numerically stable for ill-conditioned problems.
void least_squares_qr(cusolverDnHandle_t solver, float* d_A, float* d_b,
float* d_x, int m, int n) {
// cuSOLVER gels handles QR internally
int lwork;
cusolverDnSgels_bufferSize(solver, m, n, 1, d_A, m, d_b, m, d_x, n, NULL, &lwork);
float* d_work;
cudaMalloc(&d_work, lwork * sizeof(float));
int* d_info;
cudaMalloc(&d_info, sizeof(int));
// Solve min ||Ax - b|| via QR
cusolverDnSgels(solver, m, n, 1, d_A, m, d_b, m, d_x, n, d_work, lwork, d_info);
}
// Alternative: explicit QR for multiple RHS
void least_squares_qr_explicit(cusolverDnHandle_t solver, cublasHandle_t blas,
float* d_A, float* d_B, float* d_X, int m, int n, int nrhs) {
// QR factorization
float* d_tau;
cudaMalloc(&d_tau, n * sizeof(float));
cusolverDnSgeqrf(solver, m, n, d_A, m, d_tau, ...);
// Apply Q^T to B
cusolverDnSormqr(solver, CUBLAS_SIDE_LEFT, CUBLAS_OP_T,
m, nrhs, n, d_A, m, d_tau, d_B, m, ...);
// Solve R * X = Q^T * B (triangular solve)
cublasSrsm(blas, CUBLAS_SIDE_LEFT, CUBLAS_FILL_MODE_UPPER,
CUBLAS_OP_N, CUBLAS_DIAG_NON_UNIT, n, nrhs,
&one, d_A, m, d_B, m);
}| Metric | Naive | Optimized | Improvement |
|---|---|---|---|
| 10000x100 well-conditioned | 45ms (normal) | 32ms (normal+Cholesky) | 1.4x faster |
| 10000x100 ill-conditioned | Wrong answer | 85ms (QR, correct) | Correctness |
| Batch 1000 100x10 systems | 120ms (sequential) | 8ms (batched) | 15x faster |
Use normal equations when κ(A) < 1e4 and speed matters. Use QR when numerical accuracy is critical or A is ill-conditioned. For machine learning with well-scaled features, normal equations often suffice.
For ridge: solve (A^TA + λI)x = A^Tb. Add λ to diagonal of A^TA before Cholesky. For QR approach, augment A with sqrt(λ)*I and b with zeros.
Ready to optimize your CUDA code? Download RightNow AI and get real-time performance analysis for your kernels.