| 5834 | def QR(self, full=True): |
| 5835 | r""" |
| 5836 | Returns a factorization of ``self`` as a unitary matrix |
| 5837 | and an upper-triangular matrix. |
| 5838 | |
| 5839 | INPUT: |
| 5840 | |
| 5841 | - ``full`` - default: ``True`` - if ``True`` then the |
| 5842 | returned matrices have dimensions as described below. |
| 5843 | If ``False`` the ``R`` matrix has no zero rows and the |
| 5844 | columns of ``Q`` are a basis for the column space of |
| 5845 | ``self``. |
| 5846 | |
| 5847 | OUTPUT: |
| 5848 | |
| 5849 | If ``self`` is an `m\times n` matrix and ``full=True`` then this |
| 5850 | method returns a pair of matrices: `Q` is an `m\times m` unitary |
| 5851 | matrix (meaning its inverse is its conjugate-transpose) and `R` |
| 5852 | is an `m\times n` upper-triangular matrix with non-negative entries |
| 5853 | on the diagonal. For a matrix of full rank this factorization is |
| 5854 | unique (due to the restriction to positive entries on the diagonal). |
| 5855 | |
| 5856 | If ``full=False`` then `Q` has `m` rows and the columns form an |
| 5857 | orthonormal basis for the column space of ``self``. So, in particular, |
| 5858 | the conjugate-transpose of `Q` times `Q` will be an identity matrix. |
| 5859 | The matrix `R` will still be upper-triangular but will also have full |
| 5860 | rank, in particular it will lack the zero rows present in a full |
| 5861 | factorization of a rank-deficient matrix. |
| 5862 | |
| 5863 | The results obtained when ``full=True`` are cached, |
| 5864 | hence `Q` and `R` are immutable matrices in this case. |
| 5865 | |
| 5866 | .. note:: |
| 5867 | |
| 5868 | This is an exact computation, so limited to exact rings. |
| 5869 | Also the base ring needs to have a fraction field implemented |
| 5870 | in Sage and this field must contain square roots. One example |
| 5871 | is the field of algebraic numbers, ``QQbar``, as used in the |
| 5872 | examples below. If you need numerical results, convert the |
| 5873 | base ring to the field of complex double numbers, ``CDF``, |
| 5874 | which will use a faster routine that is careful about |
| 5875 | numerical subtleties. |
| 5876 | |
| 5877 | ALGORITHM: |
| 5878 | |
| 5879 | "Modified Gram-Schmidt," Algorithm 8.1 of [TREFETHEN-BAU]_. |
| 5880 | |
| 5881 | EXAMPLES: |
| 5882 | |
| 5883 | For a nonsingular matrix, the QR decomposition is unique. :: |
| 5884 | |
| 5885 | sage: A = matrix(QQbar, [[-2, 0, -4, -1, -1], |
| 5886 | ... [-2, 1, -6, -3, -1], |
| 5887 | ... [1, 1, 7, 4, 5], |
| 5888 | ... [3, 0, 8, 3, 3], |
| 5889 | ... [-1, 1, -6, -6, 5]]) |
| 5890 | sage: Q, R = A.QR() |
| 5891 | sage: Q |
| 5892 | [ -0.4588314677411235? -0.1260506983326509? 0.3812120831224489? -0.3945737113384181? -0.6874400625963308?] |
| 5893 | [ -0.4588314677411235? 0.4726901187474409? -0.05198346588033394? 0.7172941251646595? -0.2209628772631064?] |
| 5894 | [ 0.2294157338705618? 0.6617661662464172? 0.6619227988762521? -0.1808720937375480? 0.1964114464560946?] |
| 5895 | [ 0.6882472016116853? 0.1890760474989764? -0.2044682991293135? 0.0966302966543065? -0.6628886317893191?] |
| 5896 | [ -0.2294157338705618? 0.5357154679137663? -0.6099393329959182? -0.5364220314271115? 0.02455143080701182?] |
| 5897 | sage: R |
| 5898 | [ 4.358898943540674? -0.4588314677411235? 13.07669683062202? 6.194224814505168? 2.982404540317303?] |
| 5899 | [ 0 1.670171752907625? 0.5987408170800917? -1.292019657909672? 6.207996892883057?] |
| 5900 | [ 0 0 5.444401659866974? 5.468660610611130? -0.6827161852283857?] |
| 5901 | [ 0 0 0 1.027626039419836? -3.619300149686620?] |
| 5902 | [ 0 0 0 0 0.02455143080701182?] |
| 5903 | sage: Q.conjugate_transpose()*Q |
| 5904 | [1.000000000000000? 0.?e-37 0.?e-36 0.?e-34 0.?e-31] |
| 5905 | [ 0.?e-37 1.000000000000000? 0.?e-36 0.?e-34 0.?e-31] |
| 5906 | [ 0.?e-36 0.?e-36 1.000000000000000? 0.?e-34 0.?e-31] |
| 5907 | [ 0.?e-34 0.?e-34 0.?e-34 1.000000000000000? 0.?e-31] |
| 5908 | [ 0.?e-31 0.?e-31 0.?e-31 0.?e-31 1.000000000000000?] |
| 5909 | sage: Q*R == A |
| 5910 | True |
| 5911 | |
| 5912 | |
| 5913 | An example with complex numbers in ``QQbar``, the field of algebraic |
| 5914 | numbers. :: |
| 5915 | |
| 5916 | sage: A = matrix(QQbar, [[-8, 4*I + 1, -I + 2, 2*I + 1], |
| 5917 | ... [1, -2*I - 1, -I + 3, -I + 1], |
| 5918 | ... [I + 7, 2*I + 1, -2*I + 7, -I + 1], |
| 5919 | ... [I + 2, 0, I + 12, -1]]) |
| 5920 | sage: Q, R = A.QR() |
| 5921 | sage: Q |
| 5922 | [ -0.7302967433402215? 0.2070566455055649? + 0.5383472783144687?*I 0.2463049809998642? - 0.0764456358723292?*I 0.2381617683194332? - 0.10365960327796942?*I] |
| 5923 | [ 0.0912870929175277? -0.2070566455055649? - 0.3778783780476559?*I 0.3786559533863033? - 0.1952221495524667?*I 0.7012444502144683? - 0.3643711650986595?*I] |
| 5924 | [ 0.6390096504226938? + 0.0912870929175277?*I 0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I 0.3140171085506764? - 0.0825191718705412?*I] |
| 5925 | [ 0.1825741858350554? + 0.0912870929175277?*I -0.03623491296347385? + 0.0724698259269477?*I 0.8632284069415110? + 0.06322839976356195?*I -0.4499694867611521? - 0.01161191812089175?*I] |
| 5926 | sage: R |
| 5927 | [ 10.95445115010333? 0.?e-37 - 1.917028951268082?*I 5.385938482134133? - 2.190890230020665?*I -0.2738612787525831? - 2.190890230020665?*I] |
| 5928 | [ 0 4.829596256417300? -0.869637911123373? - 5.864879483945125?*I 0.993871898426712? - 0.3054085521207082?*I] |
| 5929 | [ 0 0 12.00160760935814? -0.2709533402297274? + 0.4420629644486323?*I] |
| 5930 | [ 0 0 0 1.942963944258992? + 0.?e-33*I] |
| 5931 | sage: Q.conjugate_transpose()*Q |
| 5932 | [1.000000000000000? + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-33 + 0.?e-33*I 0.?e-33 + 0.?e-33*I] |
| 5933 | [ 0.?e-36 + 0.?e-36*I 1.000000000000000? + 0.?e-36*I 0.?e-33 + 0.?e-33*I 0.?e-33 + 0.?e-33*I] |
| 5934 | [ 0.?e-33 + 0.?e-33*I 0.?e-33 + 0.?e-33*I 1.000000000000000? + 0.?e-34*I 0.?e-33 + 0.?e-33*I] |
| 5935 | [ 0.?e-33 + 0.?e-33*I 0.?e-33 + 0.?e-33*I 0.?e-33 + 0.?e-33*I 1.000000000000000? + 0.?e-33*I] |
| 5936 | sage: Q*R - A |
| 5937 | [ 0 0.?e-36 + 0.?e-36*I 0.?e-32 + 0.?e-32*I 0.?e-33 + 0.?e-33*I] |
| 5938 | [ 0 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-33 + 0.?e-33*I] |
| 5939 | [0.?e-36 + 0.?e-36*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-33 + 0.?e-33*I] |
| 5940 | [0.?e-37 + 0.?e-37*I 0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-33 + 0.?e-33*I] |
| 5941 | |
| 5942 | A rank-deficient rectangular matrix, with both values of the ``full`` keyword. :: |
| 5943 | |
| 5944 | sage: A = matrix(QQbar, [[2, -3, 3], |
| 5945 | ... [-1, 1, -1], |
| 5946 | ... [-1, 3, -3], |
| 5947 | ... [-5, 1, -1]]) |
| 5948 | sage: Q, R = A.QR() |
| 5949 | sage: Q |
| 5950 | [ 0.3592106040535498? -0.5693261797050169? 0.7239227659930268? 0.1509015305256380?] |
| 5951 | [ -0.1796053020267749? 0.1445907757980996? 0 0.9730546968377341?] |
| 5952 | [ -0.1796053020267749? 0.7048800320157352? 0.672213996993525? -0.1378927778941174?] |
| 5953 | [ -0.8980265101338745? -0.3976246334447737? 0.1551263069985058? -0.10667177157846818?] |
| 5954 | sage: R |
| 5955 | [ 5.567764362830022? -2.694079530401624? 2.694079530401624?] |
| 5956 | [ 0 3.569584777515583? -3.569584777515583?] |
| 5957 | [ 0 0 0] |
| 5958 | [ 0 0 0] |
| 5959 | sage: Q.conjugate_transpose()*Q |
| 5960 | [ 1 0.?e-37 0.?e-37 0.?e-38] |
| 5961 | [ 0.?e-37 1 0.?e-37 0.?e-38] |
| 5962 | [ 0.?e-37 0.?e-37 1 0.?e-38] |
| 5963 | [ 0.?e-38 0.?e-37 0.?e-38 1.000000000000000?] |
| 5964 | |
| 5965 | sage: Q, R = A.QR(full=False) |
| 5966 | sage: Q |
| 5967 | [ 0.3592106040535498? -0.5693261797050169?] |
| 5968 | [-0.1796053020267749? 0.1445907757980996?] |
| 5969 | [-0.1796053020267749? 0.7048800320157352?] |
| 5970 | [-0.8980265101338745? -0.3976246334447737?] |
| 5971 | sage: R |
| 5972 | [ 5.567764362830022? -2.694079530401624? 2.694079530401624?] |
| 5973 | [ 0 3.569584777515583? -3.569584777515583?] |
| 5974 | sage: Q.conjugate_transpose()*Q |
| 5975 | [ 1 0.?e-37] |
| 5976 | [0.?e-37 1] |
| 5977 | |
| 5978 | Another rank-deficient rectangular matrix, with complex entries, |
| 5979 | as a reduced decomposition. :: |
| 5980 | |
| 5981 | sage: A = matrix(QQbar, [[-3*I - 3, I - 3, -12*I + 1, -2], |
| 5982 | ... [-I - 1, -2, 5*I - 1, -I - 2], |
| 5983 | ... [-4*I - 4, I - 5, -7*I, -I - 4]]) |
| 5984 | sage: Q, R = A.QR(full=False) |
| 5985 | sage: Q |
| 5986 | [ -0.4160251471689219? - 0.4160251471689219?*I 0.5370861555295747? + 0.1790287185098583?*I] |
| 5987 | [ -0.1386750490563073? - 0.1386750490563073?*I -0.7519206177414046? - 0.2506402059138015?*I] |
| 5988 | [ -0.5547001962252291? - 0.5547001962252291?*I -0.2148344622118299? - 0.07161148740394329?*I] |
| 5989 | sage: R |
| 5990 | [ 7.211102550927979? 3.328201177351375? - 5.269651864139676?*I 7.904477796209515? + 8.45917799243475?*I 4.021576422632911? - 2.634825932069838?*I] |
| 5991 | [ 0 1.074172311059150? -1.611258466588724? - 9.13046464400277?*I 1.611258466588724? + 0.5370861555295747?*I] |
| 5992 | sage: Q.conjugate_transpose()*Q |
| 5993 | [1.000000000000000? + 0.?e-37*I 0.?e-36 + 0.?e-36*I] |
| 5994 | [ 0.?e-36 + 0.?e-36*I 1.000000000000000? + 0.?e-36*I] |
| 5995 | sage: Q*R-A |
| 5996 | [0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I] |
| 5997 | [0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-34*I 0.?e-35 + 0.?e-35*I] |
| 5998 | [0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-36*I] |
| 5999 | |
| 6000 | Results of full decompositions are cached and thus returned |
| 6001 | immutable. :: |
| 6002 | |
| 6003 | sage: A = random_matrix(QQbar, 2, 2) |
| 6004 | sage: Q, R = A.QR() |
| 6005 | sage: Q.is_mutable() |
| 6006 | False |
| 6007 | sage: R.is_mutable() |
| 6008 | False |
| 6009 | |
| 6010 | Trivial cases return trivial results of the correct size, |
| 6011 | and we check `Q` itself in one case. :: |
| 6012 | |
| 6013 | sage: A = zero_matrix(QQbar, 0, 10) |
| 6014 | sage: Q, R = A.QR() |
| 6015 | sage: Q.nrows(), Q.ncols() |
| 6016 | (0, 0) |
| 6017 | sage: R.nrows(), R.ncols() |
| 6018 | (0, 10) |
| 6019 | sage: A = zero_matrix(QQbar, 3, 0) |
| 6020 | sage: Q, R = A.QR() |
| 6021 | sage: Q.nrows(), Q.ncols() |
| 6022 | (3, 3) |
| 6023 | sage: R.nrows(), R.ncols() |
| 6024 | (3, 0) |
| 6025 | sage: Q |
| 6026 | [1 0 0] |
| 6027 | [0 1 0] |
| 6028 | [0 0 1] |
| 6029 | |
| 6030 | TESTS: |
| 6031 | |
| 6032 | Inexact rings are caught and ``CDF`` suggested. :: |
| 6033 | |
| 6034 | sage: A = matrix(RealField(100), 2, range(4)) |
| 6035 | sage: A.QR() |
| 6036 | Traceback (most recent call last): |
| 6037 | ... |
| 6038 | NotImplementedError: QR decomposition is implemented over exact rings, try CDF for numerical results, not Real Field with 100 bits of precision |
| 6039 | |
| 6040 | Without a fraction field, we cannot hope to run the algorithm. :: |
| 6041 | |
| 6042 | sage: A = matrix(Integers(6), 2, range(4)) |
| 6043 | sage: A.QR() |
| 6044 | Traceback (most recent call last): |
| 6045 | ... |
| 6046 | ValueError: QR decomposition needs a fraction field of Ring of integers modulo 6 |
| 6047 | |
| 6048 | The biggest obstacle is making unit vectors, thus requiring square |
| 6049 | roots, though some small cases pass through. :: |
| 6050 | |
| 6051 | sage: A = matrix(ZZ, 3, range(9)) |
| 6052 | sage: A.QR() |
| 6053 | Traceback (most recent call last): |
| 6054 | ... |
| 6055 | TypeError: QR decomposition unable to compute square roots in Rational Field |
| 6056 | |
| 6057 | sage: A = matrix(ZZ, 2, range(4)) |
| 6058 | sage: Q, R = A.QR() |
| 6059 | sage: Q |
| 6060 | [0 1] |
| 6061 | [1 0] |
| 6062 | sage: R |
| 6063 | [2 3] |
| 6064 | [0 1] |
| 6065 | |
| 6066 | REFERENCES: |
| 6067 | |
| 6068 | .. [TREFETHEN-BAU] Trefethen, Lloyd N., Bau, David, III |
| 6069 | "Numerical Linear Algebra" |
| 6070 | SIAM, Philadelphia, 1997. |
| 6071 | |
| 6072 | AUTHOR: |
| 6073 | |
| 6074 | - Rob Beezer (2011-02-17) |
| 6075 | """ |
| 6076 | from sage.modules.free_module_element import zero_vector |
| 6077 | from sage.matrix.constructor import zero_matrix, matrix |
| 6078 | from sage.functions.other import sqrt |
| 6079 | |
| 6080 | if full: |
| 6081 | QR = self.fetch('QR_factors') |
| 6082 | if QR is not None: |
| 6083 | return QR |
| 6084 | R = self.base_ring() |
| 6085 | if not R.is_exact(): |
| 6086 | raise NotImplementedError('QR decomposition is implemented over exact rings, try CDF for numerical results, not %s' % R) |
| 6087 | try: |
| 6088 | F = R.fraction_field() |
| 6089 | except: |
| 6090 | raise ValueError("QR decomposition needs a fraction field of %s" % R) |
| 6091 | m = self.nrows() |
| 6092 | n = self.ncols() |
| 6093 | |
| 6094 | R = zero_matrix(F, m, n) |
| 6095 | V = self.columns(copy=True) |
| 6096 | Q = [] |
| 6097 | row = 0 # row of R being filled |
| 6098 | for i in range(n): |
| 6099 | v = V[i] |
| 6100 | hip = v.hermitian_inner_product(v) |
| 6101 | if hip != 0: |
| 6102 | try: |
| 6103 | scale = sqrt(hip) |
| 6104 | q = (1/scale)*v |
| 6105 | Q.append(q) |
| 6106 | R[row,i] = scale |
| 6107 | for j in range(i+1, n): |
| 6108 | R[row,j] = q.hermitian_inner_product(V[j]) |
| 6109 | V[j] = V[j] - R[row,j]*q |
| 6110 | row = row + 1 |
| 6111 | except TypeError: |
| 6112 | raise TypeError('QR decomposition unable to compute square roots in %s' % F) |
| 6113 | # complete to full orthonormal basis, or reduce to truncated R |
| 6114 | if full: |
| 6115 | Qt = matrix(Q) # as rows here |
| 6116 | if Qt.nrows() == 0: |
| 6117 | Qt = zero_matrix(F, 0, m) |
| 6118 | orthogonal = Qt.right_kernel().basis_matrix().transpose() |
| 6119 | Qperp, _ = orthogonal.QR(full=False) |
| 6120 | Q = Q + Qperp.columns() |
| 6121 | else: |
| 6122 | R = R[0:len(Q), 0:n] |
| 6123 | Q = matrix(Q).transpose() |
| 6124 | # Adjust rows of Q if empty |
| 6125 | if Q.ncols() == 0: |
| 6126 | Q = zero_matrix(F, m, 0) |
| 6127 | QR = (Q, R) |
| 6128 | if full: |
| 6129 | Q.set_immutable() |
| 6130 | R.set_immutable() |
| 6131 | self.cache('QR_factors', QR) |
| 6132 | return QR |
| 6133 | |