只是为了好玩,我添加了此代码的犰狳版本并对其进行了基准测试
犰狳代码:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
//[[Rcpp::export]]
List grahm_schimdtCpp(arma::mat A) {
int n = A.n_cols;
int m = A.n_rows;
arma::mat Q(m, n);
Q.fill(0);
arma::mat R(n, n);
R.fill(0);
for (int j = 0; j < n; j++) {
arma::vec v = A.col(j);
if (j > 0) {
for(int i = 0; i < j; i++) {
R(i, j) = arma::as_scalar(Q.col(i).t() * A.col(j));
v = v - R(i, j) * Q.col(i);
}
}
R(j, j) = arma::norm(v, 2);
Q.col(j) = v / R(j, j);
}
return List::create(_["Q"] = Q,
_["R"] = R
);
}
R代码未优化(直接基于算法)
grahm_schimdtR <- function(A) {
m <- nrow(A)
n <- ncol(A)
Q <- matrix(0, nrow = m, ncol = n)
R <- matrix(0, nrow = n, ncol = n)
for (j in 1:n) {
v <- A[ , j, drop = FALSE]
if (j > 1) {
for(i in 1:(j-1)) {
R[i, j] <- t(Q[,i,drop = FALSE]) %*% A[ , j, drop = FALSE]
v <- v - R[i, j] * Q[ ,i]
}
}
R[j, j] = norm(v, type = "2")
Q[ ,j] = v / R[j, j]
}
list("Q" = Q, "R" = R)
}
R 中的本机 QR 分解
qrNative <- function(A) {
qrdec <- qr(A)
list(Q = qr.R(qrdec), R = qr.Q(qrdec))
}
我们将使用与原始文档中相同的矩阵来测试它(上面帖子中的链接)
A <- matrix(c(4, 3, -2, 1), ncol = 2)
all.equal(grahm_schimdtR(A)$Q %*% grahm_schimdtR(A)$R, A)
## [1] TRUE
all.equal(grahm_schimdtCpp(A)$Q %*% grahm_schimdtCpp(A)$R, A)
## [1] TRUE
all.equal(qrNative(A)$Q %*% qrNative(A)$R, A)
## [1] TRUE
现在让我们对其进行基准测试
require(rbenchmark)
set.seed(123)
A <- matrix(rnorm(10000), 100, 100)
benchmark(qrNative(A),
grahm_schimdtR(A),
grahm_schimdtCpp(A),
order = "elapsed")
## test replications elapsed relative user.self
## 3 grahm_schimdtCpp(A) 100 0.272 1.000 0.272
## 1 qrNative(A) 100 1.013 3.724 1.144
## 2 grahm_schimdtR(A) 100 84.279 309.849 95.042
## sys.self user.child sys.child
## 3 0.000 0 0
## 1 0.872 0 0
## 2 72.577 0 0
我真的很喜欢将代码移植到 Rcpp 中是多么容易......