# Introduction

This package implements in R the LAPACK routine DGEQRF to obtain the QR factorization without pivoting of a real matrix.

This package aims to solve one of the issues of the qr() function of base R, it returns the solution obtained using QR factorization with pivoting. Therefore, the product Q*R is not equal to the factorized matrix.

# Further information

For further explanation on the DGEQRF routine, read the following document.

# Example on how to use QR

First, we need to define the matrix we want to factorize. This matrix can be of any m x n dimensions and the entries must be real numbers.

# Let's sample a random square-matrix
A<-matrix(runif(121,min = -100, max = 100), 11, 11)

Once, we have the matrix, it is time to use the QR function:

QRres<-QR(A)
QRres
## $qr ## [,1] [,2] [,3] [,4] [,5] ## [1,] -187.38216511 37.48936604 -25.79571476 -18.36028428 -18.89303938 ## [2,] 0.02299492 -161.31509872 24.98486513 -99.92201658 60.16808487 ## [3,] 0.13447582 -0.03904225 -200.61442852 -30.44961762 -49.71367247 ## [4,] 0.38768007 -0.16819982 -0.17136347 182.64577465 54.13845072 ## [5,] -0.05374256 0.29586923 -0.12536500 -0.07338048 -132.78595306 ## [6,] 0.25717315 -0.06575008 0.31627776 -0.37566077 -0.07332238 ## [7,] -0.36427200 0.47219424 0.26683708 0.57437987 -0.11992795 ## [8,] -0.08603554 -0.37617664 -0.43744119 -0.05251349 0.60690638 ## [9,] 0.28734800 -0.25288232 0.12916952 0.33455781 0.45183990 ## [10,] -0.24840642 0.18694424 -0.27298415 0.06737860 0.28449435 ## [11,] 0.30101551 0.14087384 -0.08674715 0.52404056 -0.08933874 ## [,6] [,7] [,8] [,9] [,10] ## [1,] 77.2843566 -91.95513661 29.0959221 -39.58683992 87.6646902 ## [2,] -9.0884507 -102.58369054 -119.0385751 -88.45746100 -27.9144422 ## [3,] 34.3805674 -46.32264701 -42.3637496 2.55510600 -28.8193252 ## [4,] -84.4097535 -65.15015262 79.6735475 22.55690830 -33.7666197 ## [5,] -37.8882086 17.93850552 24.2436619 92.50465285 -75.1457953 ## [6,] -154.1964105 59.33751257 10.3214249 87.05428590 20.7839009 ## [7,] -0.2383296 124.52215829 -63.9796803 -47.37083659 -24.0745709 ## [8,] -0.1221353 0.30400182 -110.4494378 104.51983943 -5.3854435 ## [9,] 0.3798313 0.34333956 -0.1456755 53.52076381 -57.6825134 ## [10,] -0.1532163 0.20915181 0.3471539 0.21911489 116.2258071 ## [11,] -0.2140989 0.02157955 -0.5449169 -0.07852387 -0.6273862 ## [,11] ## [1,] -157.51426 ## [2,] -56.65363 ## [3,] -24.60077 ## [4,] -43.61615 ## [5,] 53.69243 ## [6,] 42.85188 ## [7,] 53.25999 ## [8,] 19.23805 ## [9,] -40.42576 ## [10,] -65.40451 ## [11,] -30.75442 ## ##$qraux
##  [1] 1.239990 1.246186 1.327668 1.069368 1.189646 1.556053 1.594249 1.390171
##  [9] 1.897214 1.435118 0.000000
##
## $Q ## [,1] [,2] [,3] [,4] [,5] [,6] ## [1,] -0.23999033 -0.48685847 0.2168492 0.07945991 -0.04444925 0.43725058 ## [2,] -0.02851347 -0.25738157 0.2629875 0.37364373 -0.37695819 -0.45783732 ## [3,] -0.16674872 -0.01681677 -0.3085802 0.04567595 -0.53942519 0.52901837 ## [4,] -0.48071954 0.02086298 0.2681862 -0.10958652 0.18904712 -0.25870222 ## [5,] 0.06664026 -0.34254316 0.2311237 0.17800301 -0.23337398 -0.04746951 ## [6,] -0.31889222 -0.04327008 -0.3811077 0.41336571 -0.09129722 -0.11880020 ## [7,] 0.45169375 -0.41109309 -0.3114366 -0.45438872 -0.13669904 -0.15218332 ## [8,] 0.10668324 0.51067331 0.4650660 -0.11220525 -0.33956815 0.23801999 ## [9,] -0.35630874 0.17524067 -0.1744270 -0.42256392 -0.50960267 -0.36929535 ## [10,] 0.30802156 -0.11202858 0.3567975 -0.03579657 -0.24474623 -0.07026738 ## [11,] -0.37325633 -0.32210701 0.2167920 -0.48838910 0.11323025 0.12211864 ## [,7] [,8] [,9] [,10] [,11] ## [1,] -0.48867847 -0.13839733 -0.04336296 -0.33911311 -0.28988567 ## [2,] -0.11488180 -0.05495472 -0.13858244 0.52393944 -0.25004334 ## [3,] 0.37999924 -0.29634025 -0.10457290 0.23645485 0.06856165 ## [4,] 0.30290541 -0.62979168 0.28872854 -0.09751158 -0.02475552 ## [5,] 0.57107690 0.41989794 0.27962706 -0.39267998 -0.03963951 ## [6,] -0.31516023 0.09537911 0.50388757 -0.02677351 0.44139532 ## [7,] -0.09251755 -0.24144773 0.43676066 0.07746477 -0.13434958 ## [8,] -0.20859039 0.07130070 0.50468082 0.11151895 -0.10612266 ## [9,] -0.14737981 0.15904164 -0.22923502 -0.35750602 -0.09950634 ## [10,] -0.11167987 -0.28165794 -0.23102452 -0.20267876 0.71616595 ## [11,] 0.01126817 0.37410341 0.03372389 0.45036592 0.31563106 ## ##$R
##            [,1]       [,2]       [,3]      [,4]       [,5]        [,6]
##  [1,] -187.3822   37.48937  -25.79571 -18.36028  -18.89304   77.284357
##  [2,]    0.0000 -161.31510   24.98487 -99.92202   60.16808   -9.088451
##  [3,]    0.0000    0.00000 -200.61443 -30.44962  -49.71367   34.380567
##  [4,]    0.0000    0.00000    0.00000 182.64577   54.13845  -84.409754
##  [5,]    0.0000    0.00000    0.00000   0.00000 -132.78595  -37.888209
##  [6,]    0.0000    0.00000    0.00000   0.00000    0.00000 -154.196410
##  [7,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
##  [8,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
##  [9,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
## [10,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
## [11,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
##             [,7]       [,8]       [,9]      [,10]      [,11]
##  [1,]  -91.95514   29.09592 -39.586840  87.664690 -157.51426
##  [2,] -102.58369 -119.03858 -88.457461 -27.914442  -56.65363
##  [3,]  -46.32265  -42.36375   2.555106 -28.819325  -24.60077
##  [4,]  -65.15015   79.67355  22.556908 -33.766620  -43.61615
##  [5,]   17.93851   24.24366  92.504653 -75.145795   53.69243
##  [6,]   59.33751   10.32142  87.054286  20.783901   42.85188
##  [7,]  124.52216  -63.97968 -47.370837 -24.074571   53.25999
##  [8,]    0.00000 -110.44944 104.519839  -5.385443   19.23805
##  [9,]    0.00000    0.00000  53.520764 -57.682513  -40.42576
## [10,]    0.00000    0.00000   0.000000 116.225807  -65.40451
## [11,]    0.00000    0.00000   0.000000   0.000000  -30.75442

We can observe that the output is a list with 4 elements:

• qr: This element is a matrix returned by the DGEQRF routine with the same dimensions as A. The upper triangle contains the R of the decomposition and the lower triangle contains information on the Q of the decomposition (stored in compact form).

• qraux: This element is a vector returned by the DGEQRF routine of length ncol(A) which contains additional information on Q.

• Q: This element is an orthogonal matrix such that Q*R is equal to A.

• R: This element is an upper triangular matrix such that Q*R is the input matrix.

Letâ€™s print the Q and R matrices obtained using QR().

QRres$Q ## [,1] [,2] [,3] [,4] [,5] [,6] ## [1,] -0.23999033 -0.48685847 0.2168492 0.07945991 -0.04444925 0.43725058 ## [2,] -0.02851347 -0.25738157 0.2629875 0.37364373 -0.37695819 -0.45783732 ## [3,] -0.16674872 -0.01681677 -0.3085802 0.04567595 -0.53942519 0.52901837 ## [4,] -0.48071954 0.02086298 0.2681862 -0.10958652 0.18904712 -0.25870222 ## [5,] 0.06664026 -0.34254316 0.2311237 0.17800301 -0.23337398 -0.04746951 ## [6,] -0.31889222 -0.04327008 -0.3811077 0.41336571 -0.09129722 -0.11880020 ## [7,] 0.45169375 -0.41109309 -0.3114366 -0.45438872 -0.13669904 -0.15218332 ## [8,] 0.10668324 0.51067331 0.4650660 -0.11220525 -0.33956815 0.23801999 ## [9,] -0.35630874 0.17524067 -0.1744270 -0.42256392 -0.50960267 -0.36929535 ## [10,] 0.30802156 -0.11202858 0.3567975 -0.03579657 -0.24474623 -0.07026738 ## [11,] -0.37325633 -0.32210701 0.2167920 -0.48838910 0.11323025 0.12211864 ## [,7] [,8] [,9] [,10] [,11] ## [1,] -0.48867847 -0.13839733 -0.04336296 -0.33911311 -0.28988567 ## [2,] -0.11488180 -0.05495472 -0.13858244 0.52393944 -0.25004334 ## [3,] 0.37999924 -0.29634025 -0.10457290 0.23645485 0.06856165 ## [4,] 0.30290541 -0.62979168 0.28872854 -0.09751158 -0.02475552 ## [5,] 0.57107690 0.41989794 0.27962706 -0.39267998 -0.03963951 ## [6,] -0.31516023 0.09537911 0.50388757 -0.02677351 0.44139532 ## [7,] -0.09251755 -0.24144773 0.43676066 0.07746477 -0.13434958 ## [8,] -0.20859039 0.07130070 0.50468082 0.11151895 -0.10612266 ## [9,] -0.14737981 0.15904164 -0.22923502 -0.35750602 -0.09950634 ## [10,] -0.11167987 -0.28165794 -0.23102452 -0.20267876 0.71616595 ## [11,] 0.01126817 0.37410341 0.03372389 0.45036592 0.31563106 QRres$R
##            [,1]       [,2]       [,3]      [,4]       [,5]        [,6]
##  [1,] -187.3822   37.48937  -25.79571 -18.36028  -18.89304   77.284357
##  [2,]    0.0000 -161.31510   24.98487 -99.92202   60.16808   -9.088451
##  [3,]    0.0000    0.00000 -200.61443 -30.44962  -49.71367   34.380567
##  [4,]    0.0000    0.00000    0.00000 182.64577   54.13845  -84.409754
##  [5,]    0.0000    0.00000    0.00000   0.00000 -132.78595  -37.888209
##  [6,]    0.0000    0.00000    0.00000   0.00000    0.00000 -154.196410
##  [7,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
##  [8,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
##  [9,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
## [10,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
## [11,]    0.0000    0.00000    0.00000   0.00000    0.00000    0.000000
##             [,7]       [,8]       [,9]      [,10]      [,11]
##  [1,]  -91.95514   29.09592 -39.586840  87.664690 -157.51426
##  [2,] -102.58369 -119.03858 -88.457461 -27.914442  -56.65363
##  [3,]  -46.32265  -42.36375   2.555106 -28.819325  -24.60077
##  [4,]  -65.15015   79.67355  22.556908 -33.766620  -43.61615
##  [5,]   17.93851   24.24366  92.504653 -75.145795   53.69243
##  [6,]   59.33751   10.32142  87.054286  20.783901   42.85188
##  [7,]  124.52216  -63.97968 -47.370837 -24.074571   53.25999
##  [8,]    0.00000 -110.44944 104.519839  -5.385443   19.23805
##  [9,]    0.00000    0.00000  53.520764 -57.682513  -40.42576
## [10,]    0.00000    0.00000   0.000000 116.225807  -65.40451
## [11,]    0.00000    0.00000   0.000000   0.000000  -30.75442

We can check if Q is orthogonal by multiplying Q and $$Q^T$$.

all.equal(QRres$Q%*%t(QRres$Q),diag(11))
## [1] TRUE

Also, we can easily check that R is upper triangular.

all.equal(QRres$R[lower.tri(QRres$R)],rep(0,11*10/2))
## [1] TRUE

Finally, we can test if the factorization is correct by multiplying Q and R to obtain A.

all.equal(QRres$Q%*%QRres$R,A)
## [1] TRUE

# Comparison between QR() and qr() from base R.

In the help file of the functions to reconstruct the Q and R matrices from base R we find the following example:

# example of pivoting
x <- cbind(int = 1,
b1 = rep(1:0, each = 3), b2 = rep(0:1, each = 3),
c1 = rep(c(1,0,0), 2), c2 = rep(c(0,1,0), 2), c3 = rep(c(0,0,1),2))
x # is singular, columns "b2" and "c3" are "extra"
##      int b1 b2 c1 c2 c3
## [1,]   1  1  0  1  0  0
## [2,]   1  1  0  0  1  0
## [3,]   1  1  0  0  0  1
## [4,]   1  0  1  1  0  0
## [5,]   1  0  1  0  1  0
## [6,]   1  0  1  0  0  1
a <- qr(x)

If we multiply the obtained Q and R, we do not obtain x:

all.equal(qr.Q(a)%*%qr.R(a),x)
## [1] "Attributes: < Component \"dimnames\": Component 2: 3 string mismatches >"
## [2] "Mean relative difference: 1"

This is caused because the solution was obtained using column pivoting.

If we use the function QR(), the product of Q and R is equal to x:

all.equal(QR(x)$Q%*%QR(x)$R,x)
## [1] "Attributes: < Length mismatch: comparison on first 1 components >"

Note: The unpivoted solution is NOT equal to the pivoted solution after a transformation that undoes the column pivoting.