基于Givens的QR分解：矩形矩阵实现优化及适配工具寻求

I get it—having a working Givens QR function that handles rectangular matrices is great, but that 134x speed gap compared to pracma::givens() is frustrating, especially since pracma only supports square matrices. Let's break down how to speed up your code, or find robust alternatives that handle rectangular inputs efficiently.

First: Why Your Current Code Is Slow #

Your implementation uses sparse matrices and nested sapply/vapply loops with global assignments (<<-), which are major bottlenecks in R:

• Sparse matrix operations have unnecessary overhead for small, dense modifications (Givens rotations only touch 4 elements, so sparse matrices are overkill here).
• Global assignments (<<-) force R to look up and modify variables outside the loop scope, which is slow in repeated iterations.
• R's interpreted loops (even with vapply) can't compete with the compiled C/Fortran code underpinning pracma::givens().

Here's your original implementation for reference:

Givens.fn<-function(V) { 
  tol<-1.0*10^-14 
  m<-dim(V)[1] 
  n<-dim(V)[2] 
  spId<-bandSparse(m,m,0,list(rep(1, m+1))) # computed just once 
  Q<-spId 
  R<-V 
  sapply(1:n, function(j){ 
    if (j<m) { 
      vapply((j+1):m, function(i){ 
        # vectorized internal loop 
        if(abs(R[i,j])>tol) { 
          G<-spId 
          x<-R[j,j] 
          y<-R[i,j] 
          norm<-sqrt(x^2+y^2) 
          c<-x/norm 
          s<-y/norm 
          G[j,j]=c 
          G[i,i]=c 
          G[j,i]=s 
          G[i,j]=-s 
          Q<<-G%*%Q 
          R<<-G%*%R 
        } 
        return(1) # --> saves 15% execution time! 
      },FUN.VALUE=1.0) 
    } 
  }) 
  Q<-t(Q) 
  return(list(Q,R)) 
}

And your benchmark results (super helpful for context):

m<-20; n<-20 
set.seed(1) 
X <- as.matrix(replicate(n, runif(m))) 
library(rbenchmark) 
library(pracma) 
benchmark(Givens.fn(X), givens(X), order = "elapsed", replications = 10)

test replications elapsed relative user.self sys.self user.child sys.child
2 givens(X)          10    0.08    1.000      0.08        0         NA        NA
1 Givens.fn(X)       10   10.77  134.625     10.69      0.08         NA        NA

Optimization 1: Rewrite to Avoid Sparse Matrices & Global Assignments #

The biggest win here is to ditch sparse matrices and update Q and R directly instead of creating a full G matrix each time. Givens rotations only affect 2 rows (for R) and 2 columns (for Q), so we can compute those updates in-place.

Here's a simplified, faster version using dense matrices and local operations:

givens_qr_fast <- function(V) {
  tol <- 1e-14
  m <- nrow(V)
  n <- ncol(V)
  Q <- diag(m)
  R <- V
  
  for (j in 1:n) {
    if (j >= m) break
    for (i in (j+1):m) {
      y <- R[i, j]
      if (abs(y) <= tol) next
      x <- R[j, j]
      norm <- sqrt(x^2 + y^2)
      c <- x / norm
      s <- y / norm
      
      # Update R: apply G to rows j and i
      R_j <- R[j, ]
      R_i <- R[i, ]
      R[j, ] <- c * R_j + s * R_i
      R[i, ] <- -s * R_j + c * R_i
      
      # Update Q: apply G^T to columns j and i
      Q_j <- Q[, j]
      Q_i <- Q[, i]
      Q[, j] <- c * Q_j + s * Q_i
      Q[, i] <- -s * Q_j + c * Q_i
    }
  }
  
  # Zero out tiny elements in R to match standard QR output
  R[abs(R) < tol] <- 0
  list(Q = Q, R = R)
}

This version cuts out sparse matrix overhead and global assignments, which should already give you a massive speed boost compared to your original code.

Optimization 2: Use Compiled Code (Rcpp) for Maximum Speed #

If you need even faster performance (approaching pracma speeds), rewrite the core loop in C++ using Rcpp. Compiled code eliminates R's loop overhead entirely.

Here's a basic Rcpp implementation for rectangular matrices (install the Rcpp package first):

#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;

// [[Rcpp::export]]
List givens_qr_rcpp(NumericMatrix V) {
  int m = V.nrow();
  int n = V.ncol();
  NumericMatrix Q(m, m);
  NumericMatrix R = clone(V);
  
  // Initialize Q as identity matrix
  for (int i = 0; i < m; ++i) {
    Q(i, i) = 1.0;
  }
  
  double tol = 1e-14;
  for (int j = 0; j < n; ++j) {
    if (j >= m) break;
    for (int i = j + 1; i < m; ++i) {
      double y = R(i, j);
      if (std::abs(y) <= tol) continue;
      double x = R(j, j);
      double norm = std::sqrt(x*x + y*y);
      double c = x / norm;
      double s = y / norm;
      
      // Update R rows j and i
      for (int k = j; k < n; ++k) {
        double r_j = R(j, k);
        double r_i = R(i, k);
        R(j, k) = c * r_j + s * r_i;
        R(i, k) = -s * r_j + c * r_i;
      }
      
      // Update Q columns j and i
      for (int k = 0; k < m; ++k) {
        double q_j = Q(k, j);
        double q_i = Q(k, i);
        Q(k, j) = c * q_j + s * q_i;
        Q(k, i) = -s * q_j + c * q_i;
      }
    }
  }
  
  // Zero out small elements
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) {
      if (std::abs(R(i, j)) < tol) R(i, j) = 0.0;
    }
  }
  
  return List::create(Named("Q") = Q, Named("R") = R);
}

To use this, save it as givens_qr.cpp, then run:

library(Rcpp)
sourceCpp("givens_qr.cpp")

This should get you within spitting distance of pracma's speed, while fully supporting rectangular matrices.

Alternative: Use Existing Linear Algebra Libraries #

If you don't want to roll your own, consider using RcppEigen or RcppArmadillo—these packages wrap high-performance C++ linear algebra libraries that include fast QR decomposition (while they use Householder transformations instead of Givens, they still produce valid QR results for rectangular matrices).

For example, with RcppEigen:

library(RcppEigen)

eigen_qr <- function(V) {
  Eigen::Map<Eigen::MatrixXd> mat(as.matrix(V), nrow(V), ncol(V))
  Eigen::HouseholderQR<Eigen::MatrixXd> qr(mat)
  Q <- qr.householderQ()
  R <- qr.matrixQR().triangularView<Eigen::Upper>()
  list(Q = as.matrix(Q), R = as.matrix(R))
}

This is extremely fast and requires no manual loop writing.

内容的提问来源于stack exchange，提问作者Antonio Piemontese