Gram Schmidt with R

Question

Here is a MATLAB code for performing Gram Schmidt in page 1 http://web.mit.edu/18.06/www/Essays/gramschmidtmat.pdf

I am trying for hours and hours to perform this with R since I don't have MATLAB Here is my R

f=function(x){
m=nrow(x);
n=ncol(x);
Q=matrix(0,m,n);
R=matrix(0,n,n);

for(j in 1:n){
v=x[,j,drop=FALSE];

for(i in 1:j-1){
R[i,j]=t(Q[,i,drop=FALSE])%*%x[,j,drop=FALSE];
v=v-R[i,j]%*%Q[,i,drop=FALSE]
}

R[j,j]=max(svd(v)$d);
Q[,j,,drop=FALSE]=v/R[j,j]}

return(list(Q,R))}

It keeps on saying there is errors in either:

v=v-R[i,j]%*%Q[,i,drop=FALSE] 

or

R[j,j]=max(svd(v)$d);

What is it that I am doing wrong translating MATLAB code to R???

Answer 1

Just for fun I added an Armadillo version of this code and benchmark it

Armadillo code :

#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 code not optimized (directly based on algorithm)

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)

}

Native QR decomposition in R

qrNative <- function(A) {
    qrdec <- qr(A)
    list(Q = qr.R(qrdec), R = qr.Q(qrdec))
}

We will test it with the same matrix as in original document (link in the post above)

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

Now let's benchmark it

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

I really love how easy to port code into Rcpp....

Answer 2

If you are translating code in Matlab into R, then code semantics (code logic) should remain same. For example, in your code, you are transposing Q in t(Q[,i,drop=FALSE]) as per the given Matlab code. But Q[,i,drop=FALSE] does not return the column in column vector. So, we can make it a column vector by using the statement:

matrix(Q[,i],n,1); # n is the number of rows.

There is no error in R[j,j]=max(svd(v)$d) if v is a vector (row or column).

Yes, there is an error in

v=v-R[i,j]%*%Q[,i,drop=FALSE]

because you are using a matrix multiplication. Instead you should use a normal multiplication:

v=v-R[i,j] * Q[,i,drop=FALSE]

Here R[i,j] is a number, whereas Q[,i,drop=FALSE] is a vector. So, dimension mismatch arises here.

One more thing, if j is 3 , then 1:j-1 returns [0,1,2]. So, it should be changed to 1:(j-1), which returns [1,2] for the same value for j. But there is a catch. If j is 2, then 1:(j-1) returns [1,0]. So, 0th index is undefined for a vector or a matrix. So, we can bypass 0 value by putting a conditional expression.

Here is a working code for Gram Schmidt algorithm:

A = matrix(c(4,3,-2,1),2,2)
m = nrow(A)
n = ncol(A)
Q = matrix(0,m,n)
R = matrix(0,n,n)

for(j in 1:n)
{
    v = matrix(A[,j],n,1)
    for(i in 1:(j-1))
    {
        if(i!=0)
        {
            R[i,j] = t(matrix(Q[,i],n,1))%*%matrix(A[,j],n,1)
            v = v - (R[i,j] * matrix(Q[,i],n,1))
        }
    }
    R[j,j] = svd(v)$d 
    Q[,j] = v/R[j,j]
}

If you need to wrap the code into a function, you can do so as per your convenience.