Thunder Schmidt orthonormalization (differently) using Rcpp

Question

I am trying to get a matrix consisting of ordinary normalized orthogonal vectors based on an m-on-n-matrix, where m <n using Rcpp inside R. Unfortunately, I could not figure out how to use the qr method implemented in R to do this is.

My attempt closely follows the standard Gram-Schmidt orthonormalization algorithm, which was discussed here. However, the approach presented in this thread using Rcpp and Armadillo no longer works.

Matlab's implementation of the problem is provided by Gary Kuop in his collection of computer codes :

function [Q] = grams(A) [m,~] = size(A); Q=0*A; % compute QR using Gram-Schmidt for j = 1:m v = A(j,:); for i=1:j-1 v = v - Q(i,:)*A(j,:)'*Q(i,:); end Q(j,:) = v/norm(v); end 

My current Rcpp code is as follows:

 #include <RcppArmadillo.h> // [[Rcpp::depends(RcppArmadillo)]] using namespace Rcpp; // [[Rcpp::export]] SEXP gramscpp(arma::mat A){ int m = A.n_rows; int n = A.n_cols; arma::mat Q(m,n); Q.fill(0); Q.row(0) = A.row(0)/arma::norm(A.row(0),2); for(int j = 1; j < m; j++){ arma::rowvec v = A.row(j); for(int i = 0; i < (j-1); i++){ v = v - (Q.row(i) * A.row(j).t()) * Q.row(i); } Q.row(j) = v/arma::norm(v,2); } return Rcpp::wrap(Q); } 

When modeling a random m-by-n matrix and applying an algorithm of errors, no errors arise. However, the resulting matrix, multiplied by its transposition, does not give a unit matrix, and I can not identify the error.