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.
For further explanation on the DGEQRF routine, read the following document.
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.75442We 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.75442We can check if Q is orthogonal by multiplying Q and \(Q^T\).
all.equal(QRres$Q%*%t(QRres$Q),diag(11))
## [1] TRUEAlso, we can easily check that R is upper triangular.
all.equal(QRres$R[lower.tri(QRres$R)],rep(0,11*10/2))
## [1] TRUEFinally, we can test if the factorization is correct by multiplying Q and R to obtain A.
all.equal(QRres$Q%*%QRres$R,A)
## [1] TRUEIn 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.