QR factorization without pivoting

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 QT.

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.