This package implements in R the LAPACK routine DGEQRF to obtain theQR factorization without pivoting of a real matrix.
This package aims to solve one of the issues of the qr() function ofbase R, it returns the solution obtained using QR factorization withpivoting. Therefore, the product Q*R is not equal to the factorizedmatrix.
For further explanation on the DGEQRF routine, read the followingdocument.
First, we need to define the matrix we want to factorize. This matrixcan be of any m x n dimensions and the entries must be real numbers.
# Let's sample a random square-matrixA<-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 withthe same dimensions as A. The upper triangle contains theR of the decomposition and the lower triangle containsinformation on theQ of the decomposition (stored incompact form).
qraux: This element is a vector returned by the DGEQRF routine oflength ncol(A) which contains additional information onQ.
Q: This element is an orthogonal matrix such that Q*R is equal toA.
R: This element is an upper triangular matrix such that Q*R isthe 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.31563106QRres$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 Qand 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 matricesfrom base R we find the following example:
# example of pivotingx<-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 1a<-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 columnpivoting.
If we use the function QR(), the product of Q and R is equal tox:
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 thepivoted solution after a transformation that undoes the columnpivoting.