2234 | | sage: m = matrix(RDF,3,range(0, 12)); m |
2235 | | [ 0.0 1.0 2.0 3.0] |
2236 | | [ 4.0 5.0 6.0 7.0] |
2237 | | [ 8.0 9.0 10.0 11.0] |
2238 | | sage: Q,R = m.QR() |
2239 | | sage: Q*R |
2240 | | [ 0.0 1.0 2.0 3.0] |
2241 | | [ 4.0 5.0 6.0 7.0] |
2242 | | [ 8.0 9.0 10.0 11.0] |
2243 | | |
2244 | | Note that Q is an orthogonal matrix:: |
2245 | | |
2246 | | sage: (Q*Q.transpose()).zero_at(1e-10) |
| 2243 | Calls "linalg.qr" from SciPy, which is in turn an |
| 2244 | interface to LAPACK routines. |
| 2245 | |
| 2246 | EXAMPLES: |
| 2247 | |
| 2248 | Over the reals, the inverse of `Q` is its transpose, |
| 2249 | since including a conjugate has no effect. In the real |
| 2250 | case, we say `Q` is orthogonal. :: |
| 2251 | |
| 2252 | sage: A = matrix(RDF, [[-2, 0, -4, -1, -1], |
| 2253 | ... [-2, 1, -6, -3, -1], |
| 2254 | ... [1, 1, 7, 4, 5], |
| 2255 | ... [3, 0, 8, 3, 3], |
| 2256 | ... [-1, 1, -6, -6, 5]]) |
| 2257 | sage: Q, R = A.QR() |
| 2258 | sage: Q.round(6).zero_at(10^-6) |
| 2259 | [-0.458831 0.126051 -0.381212 -0.394574 -0.68744] |
| 2260 | [-0.458831 -0.47269 0.051983 0.717294 -0.220963] |
| 2261 | [ 0.229416 -0.661766 -0.661923 -0.180872 0.196411] |
| 2262 | [ 0.688247 -0.189076 0.204468 0.09663 -0.662889] |
| 2263 | [-0.229416 -0.535715 0.609939 -0.536422 0.024551] |
| 2264 | sage: R.round(6).zero_at(10^-6) |
| 2265 | [ 4.358899 -0.458831 13.076697 6.194225 2.982405] |
| 2266 | [ 0.0 -1.670172 -0.598741 1.29202 -6.207997] |
| 2267 | [ 0.0 0.0 -5.444402 -5.468661 0.682716] |
| 2268 | [ 0.0 0.0 0.0 1.027626 -3.6193] |
| 2269 | [ 0.0 0.0 0.0 0.0 0.024551] |
| 2270 | sage: (Q*Q.transpose()).zero_at(10^-14) |
| 2271 | [1.0 0.0 0.0 0.0 0.0] |
| 2272 | [0.0 1.0 0.0 0.0 0.0] |
| 2273 | [0.0 0.0 1.0 0.0 0.0] |
| 2274 | [0.0 0.0 0.0 1.0 0.0] |
| 2275 | [0.0 0.0 0.0 0.0 1.0] |
| 2276 | sage: (Q*R - A).zero_at(10^-14) |
| 2277 | [0.0 0.0 0.0 0.0 0.0] |
| 2278 | [0.0 0.0 0.0 0.0 0.0] |
| 2279 | [0.0 0.0 0.0 0.0 0.0] |
| 2280 | [0.0 0.0 0.0 0.0 0.0] |
| 2281 | [0.0 0.0 0.0 0.0 0.0] |
| 2282 | |
| 2283 | Now over the complex numbers, demonstrating that the SciPy libraries |
| 2284 | are (properly) using the Hermitian inner product, so that `Q` is |
| 2285 | a unitary matrix (its inverse is the conjugate-transpose). :: |
| 2286 | |
| 2287 | sage: A = matrix(CDF, [[-8, 4*I + 1, -I + 2, 2*I + 1], |
| 2288 | ... [1, -2*I - 1, -I + 3, -I + 1], |
| 2289 | ... [I + 7, 2*I + 1, -2*I + 7, -I + 1], |
| 2290 | ... [I + 2, 0, I + 12, -1]]) |
| 2291 | sage: Q, R = A.QR() |
| 2292 | sage: Q.round(6).zero_at(10^-6) |
| 2293 | [ -0.730297 0.207057 + 0.538347*I -0.246305 + 0.076446*I 0.238162 - 0.10366*I] |
| 2294 | [ 0.091287 -0.207057 - 0.377878*I -0.378656 + 0.195222*I 0.701244 - 0.364371*I] |
| 2295 | [ 0.63901 + 0.091287*I 0.170822 + 0.667758*I 0.034115 - 0.040902*I 0.314017 - 0.082519*I] |
| 2296 | [ 0.182574 + 0.091287*I -0.036235 + 0.07247*I -0.863228 - 0.063228*I -0.449969 - 0.011612*I] |
| 2297 | sage: R.round(6).zero_at(10^-6) |
| 2298 | [ 10.954451 -1.917029*I 5.385938 - 2.19089*I -0.273861 - 2.19089*I] |
| 2299 | [ 0.0 4.829596 -0.869638 - 5.864879*I 0.993872 - 0.305409*I] |
| 2300 | [ 0.0 0.0 -12.001608 0.270953 - 0.442063*I] |
| 2301 | [ 0.0 0.0 0.0 1.942964] |
| 2302 | sage: (Q.conjugate().transpose()*Q).zero_at(10^-15) |
| 2303 | [1.0 0.0 0.0 0.0] |
| 2304 | [0.0 1.0 0.0 0.0] |
| 2305 | [0.0 0.0 1.0 0.0] |
| 2306 | [0.0 0.0 0.0 1.0] |
| 2307 | sage: (Q*R - A).zero_at(10^-14) |
| 2308 | [0.0 0.0 0.0 0.0] |
| 2309 | [0.0 0.0 0.0 0.0] |
| 2310 | [0.0 0.0 0.0 0.0] |
| 2311 | [0.0 0.0 0.0 0.0] |
| 2312 | |
| 2313 | An example of a rectangular matrix that is also rank-deficient. |
| 2314 | If you run this example yourself, you may see a very small, nonzero |
| 2315 | entries in the third row, in the third column, even though the exact |
| 2316 | version of the matrix has rank 2. The final two columns of ``Q`` span the left kernel of ``A`` and different platforms will compute different |
| 2317 | bases so we do not exhibit the actual matrix. :: |
| 2318 | |
| 2319 | sage: Arat = matrix(QQ, [[2, -3, 3], |
| 2320 | ... [-1, 1, -1], |
| 2321 | ... [-1, 3, -3], |
| 2322 | ... [-5, 1, -1]]) |
| 2323 | sage: Arat.rank() |
| 2324 | 2 |
| 2325 | sage: A = Arat.change_ring(CDF) |
| 2326 | sage: Q, R = A.QR() |
| 2327 | sage: R.round(6).zero_at(10^-6) |
| 2328 | [-5.567764 2.69408 -2.69408] |
| 2329 | [ 0.0 -3.569585 3.569585] |
| 2330 | [ 0.0 0.0 0.0] |
| 2331 | [ 0.0 0.0 0.0] |
| 2332 | sage: (Q.conjugate_transpose()*Q).zero_at(10^-14) |
| 2333 | [1.0 0.0 0.0 0.0] |
| 2334 | [0.0 1.0 0.0 0.0] |
| 2335 | [0.0 0.0 1.0 0.0] |
| 2336 | [0.0 0.0 0.0 1.0] |
| 2337 | |
| 2338 | Results are cached, meaning they are immutable matrices. |
| 2339 | Make a copy if you need to manipulate a result. :: |
| 2340 | |
| 2341 | sage: A = random_matrix(CDF, 2, 2) |
| 2342 | sage: Q, R = A.QR() |
| 2343 | sage: Q.is_mutable() |
| 2344 | False |
| 2345 | sage: R.is_mutable() |
| 2346 | False |
| 2347 | sage: Q[0,0] = 0 |
| 2348 | Traceback (most recent call last): |
| 2349 | ... |
| 2350 | ValueError: matrix is immutable; please change a copy instead (i.e., use copy(M) to change a copy of M). |
| 2351 | sage: Qcopy = copy(Q) |
| 2352 | sage: Qcopy[0,0] = 679 |
| 2353 | sage: Qcopy[0,0] |
| 2354 | 679.0 |
| 2355 | |
| 2356 | TESTS: |
| 2357 | |
| 2358 | Trivial cases return trivial results of the correct size, |
| 2359 | and we check `Q` itself in one case, verifying a fix for |
| 2360 | :trac:`10795`. :: |
| 2361 | |
| 2362 | sage: A = zero_matrix(RDF, 0, 10) |
| 2363 | sage: Q, R = A.QR() |
| 2364 | sage: Q.nrows(), Q.ncols() |
| 2365 | (0, 0) |
| 2366 | sage: R.nrows(), R.ncols() |
| 2367 | (0, 10) |
| 2368 | sage: A = zero_matrix(RDF, 3, 0) |
| 2369 | sage: Q, R = A.QR() |
| 2370 | sage: Q.nrows(), Q.ncols() |
| 2371 | (3, 3) |
| 2372 | sage: R.nrows(), R.ncols() |
| 2373 | (3, 0) |
| 2374 | sage: Q |