If A is a m ×n matrix whose columns are linearly independent, the QR decomposition of A is QR = A, where Q is an m ×m orthogonal matrix and R is an upper triangular m ×n matrix.
i1 : A = random(RR^5, RR^3) o1 = | .892712 .714827 .909047 | | .673395 .89189 .314897 | | .29398 .231053 .0741835 | | .632944 .461944 .808694 | | .0258884 .775187 .362835 | 5 3 o1 : Matrix RR <--- RR 53 53 |
i2 : (Q,R) = QRDecomposition A o2 = (| -.677131 .143091 .247563 |, | -1.31837 -1.22811 -1.1883 |) | -.510777 -.324257 -.666994 | | 0 -.816018 -.175579 | | -.222987 .0524494 -.347041 | | 0 0 .523234 | | -.480095 .156448 .507738 | | -.0196366 -.92041 .339994 | o2 : Sequence |
R is upper triangular, and Q is (close to) orthogonal.
i3 : R o3 = | -1.31837 -1.22811 -1.1883 | | 0 -.816018 -.175579 | | 0 0 .523234 | 3 3 o3 : Matrix RR <--- RR 53 53 |
i4 : (transpose Q) * Q o4 = | 1 -2.77556e-17 2.77556e-17 | | -2.77556e-17 1 -2.498e-16 | | 2.77556e-17 -2.498e-16 1 | 3 3 o4 : Matrix RR <--- RR 53 53 |
i5 : clean(1e-10, oo) o5 = | 1 0 0 | | 0 1 0 | | 0 0 1 | 3 3 o5 : Matrix RR <--- RR 53 53 |
i6 : R - (transpose Q) * A o6 = | 0 0 2.22045e-16 | | 0 -3.33067e-16 0 | | 8.32667e-17 -1.66533e-16 1.11022e-16 | 3 3 o6 : Matrix RR <--- RR 53 53 |
i7 : clean(1e-10, oo) o7 = | 0 0 0 | | 0 0 0 | | 0 0 0 | 3 3 o7 : Matrix RR <--- RR 53 53 |
If the input is a MutableMatrix, then so are the output matrices.
i8 : A = mutableMatrix(RR_53, 13, 5); |
i9 : fillMatrix A o9 = | .706096 .606588 .605659 .174853 .370833 | | .127435 .848005 .96518 .626892 .339222 | | .254482 .191734 .681683 .350611 .062212 | | .741046 .403215 .914199 .379495 .465736 | | .108386 .615911 .887381 .237252 .40273 | | .348931 .0147867 .169813 .116721 .164647 | | .562428 .223028 .965004 .444183 .713493 | | .246268 .388829 .0647412 .644366 .909537 | | .153346 .557119 .877846 .194945 .566034 | | .830833 .873708 .0340514 .518585 .305423 | | .538554 .7037 .507989 .987173 .732358 | | .873665 .681869 .150294 .568273 .562839 | | .415912 .276259 .656391 .184779 .629991 | o9 : MutableMatrix |
i10 : (Q,R) = QRDecomposition A o10 = (| -.373219 .00053571 -.0145186 .434826 -.0273305 |, | -1.89191 | -.0673582 -.640249 -.160526 -.10301 .310349 | | 0 | -.134511 .0235031 -.312596 -.211448 .462328 | | 0 | -.391693 .202931 -.312943 .0456648 .146375 | | 0 | -.0572895 -.453214 -.233381 .214532 -.100994 | | 0 | -.184434 .24735 -.0733514 -.0385991 .0314994 | | -.297281 .225982 -.452339 -.214843 -.131288 | | -.130169 -.153513 .179847 -.466939 -.58274 | | -.0810536 -.368715 -.2545 .236721 -.32996 | | -.439151 -.138065 .495715 .204262 .184438 | | -.284662 -.208589 .0850958 -.577067 .107535 | | -.461791 .0602065 .325632 .0156234 -.0364137 | | -.219838 .070582 -.240391 .103952 -.379431 | ----------------------------------------------------------------------- -1.62694 -1.56271 -1.38193 -1.52288 |) -1.15333 -.947084 -.695353 -.566482 | 0 -1.57949 -.0484065 -.388761 | 0 0 -.784772 -.515486 | 0 0 0 -.778323 | o10 : Sequence |
i11 : Q*R-A o11 = | -2.22045e-16 5.55112e-16 1.11022e-15 2.77556e-16 3.33067e-16 | | 0 0 2.22045e-16 0 1.66533e-16 | | 0 5.55112e-17 -1.11022e-16 -5.55112e-17 -1.11022e-16 | | 1.11022e-16 2.22045e-16 4.44089e-16 0 1.11022e-16 | | 0 0 1.11022e-16 0 0 | | 0 1.11022e-16 1.94289e-16 5.55112e-17 -2.77556e-17 | | 0 1.66533e-16 1.11022e-16 0 1.11022e-16 | | 0 1.11022e-16 1.11022e-16 1.11022e-16 1.11022e-16 | | 0 0 1.11022e-16 2.77556e-17 0 | | 0 2.22045e-16 3.33067e-16 1.11022e-16 2.22045e-16 | | 0 1.11022e-16 2.22045e-16 -1.11022e-16 0 | | -1.11022e-16 2.22045e-16 3.33067e-16 1.11022e-16 1.11022e-16 | | -5.55112e-17 1.11022e-16 2.22045e-16 0 1.11022e-16 | o11 : MutableMatrix |
i12 : clean(1e-10,oo) o12 = | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | | 0 0 0 0 0 | o12 : MutableMatrix |
This function works by calling lapack routines, and so only uses the first 53 bits of precision. Lapack also has a way of returning an encoded pair of matrices that contain enough information to reconstruct Q, R.
If the matrices are over higher precision real or complex fields, such as RR100, this extra precision is not used in the computation