c++ - R 或 C++ 中基于 Givens 的 QR 分解：矩形矩阵的性能较差（或缺乏）实现

Question

我自己在 R 中的基于Givens的 QR 分解函数（来自都灵大学线性代数系的伪代码）如下：

>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))
}

它工作正常，但是pracma在 R 中提供的givens()函数在大约 10" 的总经过时间上快 134 倍（在我的 intel i5-3570 4core，8GB 桌面上）：

>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.00         NA        NA
1 Givens.fn(X)           10   10.77  134.625     10.69     0.08         NA        NA

不幸的是，pracma版本只接受方阵，这不是我的情况。任何进一步的代码改进？任何其他基于 Givens 的 QR 函数（甚至在 C++ 中）？谢谢。

标签： c++rshared-librarieslinear-algebra