QRGS {lazy.mat}

R Documentation

QR decomposition of X matrix by the Gram Schmidt orthogonalization or
Finding non redundant columns of X

Description

QR decomposition of X matrix by the Gram Schmidt orthogonalization or
Finding non redundant columns of X

Usage

QRGS(X, orth = 0, zero = 1, epsg = 1e-09, permute = 0, print = 0, debug = 0)

Arguments

X	The input matrix of size n x p of rank(X)=r.
orth	If orth = 0 Q will contains non redundant columns of X If orth = 1 Q will be orthogonalized (t(Q)%*%Q=diagonal) using Gram-Schumidt. If orth = 2 Q will be orthonormalized using (t(Q)%*%Q=I) Gram-Schumidt. (genuine QR decomposition)
zero	= 1 (not used)
epsg	crit for zero
permute	= 1 to enable partial pivoting (not yet available)
print	= 1 to print the summary = 2 to print Q,R,P = 3 to print Q0, P0, R0 when zero=1
debug	= 1 to print intermediate result

Details

X matrix will be decomposed as X = Q %*% R where
Q is the nrow(X) x rank(X) matrix and
R is the rank(X) x ncol(X) upper triangular matrix.
Q matrix can be obtained as Q = X %*% P.

Reduction is performed from left to right.
X may contain all zero rows/coluumns.

WHen nrow(X) > ncol(X), the Q, R and P matrices can be calculated by the native qr function as

    qr=qr(X)
    Qfromqr=qr.qy(qr,diag(nrow(X)))[,1:ncol(X)][,1:qr$rank,drop=F]
    Rfromqr=t(Qfromqr) %*% X
    Pfromqr=t(Rfromqr)%*%solve(Rfromqr%*%t(Rfromqr))
    Xhat=Qfromqr%*%Rfromqr
    QasXP=X%*%Pfromqr
    

This should always give X-Xhat == 0.

Value

A list of
Q n x r matrix, where r is the rank of X matrix
R r x p upper triangle matrix
P p x r upper triangle matrix to calculate Q as
Q = X %*% P
rank rank of X: # of independent columns of X
loczero a vector consisting of the locations of zero pivot
maxadX max( abs( X - Q %*% R ) )
maxadQ max( abs( Q - X %*% P ) )
maxadorth max(abs(offdiag(t(Q)%*%Q)))
maxadorthN max(abs(diag(t(Q)%*%Q)-I))

In addition to the above, if X matrix is not of full rank, the following matrices will be returned.
Q0 n x p matrix,
R0 p x p upper triangle matrix
P0 p x p upper triangle matrix to calculate Q0 as
Q0 = X %*% P0
These matrices contain all zero rows/columns.
maxadX0 max( abs( X - Q0 %*% R0 ) )
maxadQ0 max( abs( Q0 - X %*% P0 ) )

Examples

X=cbind(demomat(4,3),V4=c(1,3,2,5))
resQRGS <- QRGS(X, print=1)
X - resQRGS$Q %*% resQRGS$R
resQRGS$Q - X %*% resQRGS$P

resQRGS1 <- QRGS(X, orth=1, print=1)
X - resQRGS1$Q %*% resQRGS1$R
resQRGS1$Q - X %*% resQRGS1$P

resQRGS2 <- QRGS(X, orth=2, print=1)
X - resQRGS2$Q %*% resQRGS2$R
resQRGS2$Q -X %*% resQRGS2$P

# non full rank design matrix
nA <- 2
nB <- 3
# XA <- design_expand( 1:nA, type=-1, froot="A" )
# XB <- design_expand( 1:nB, type=-1, froot="B" )
XA <- model.matrix(~as.factor(1:nA) - 1)
XB <- model.matrix(~as.factor(1:nB) - 1)
Xab <- cprod( XA, XB )
colnames(Xab) <- c(paste("A",1:nA,sep=""),paste("B",1:nB,sep=""))
XAB <- Reduce( rbind
         , lapply( as.data.frame( Xab[,1:nA] )
             , function(x) t(matrix(x,nrow(Xab),nB)*Xab[,(1:nB)+nA])  ) )
colnames(XAB) <- paste("AB"
        , apply(cprod(1:nA,1:nB),1,paste,collapse=""), sep="")
X <- cbind( mu=1, Xab, XAB )

# independent columns
qr <- QRGS(X, print=3)

---