QR decomposition: Difference between revisions
Content added Content deleted
Walterpachl (talk | contribs) m (two missing 'n's inserted :-)) |
(Python solution) |
||
Line 1,309: | Line 1,309: | ||
{{works with|PARI/GP|2.6.0 and above}} |
{{works with|PARI/GP|2.6.0 and above}} |
||
<lang parigp>matqr(M)</lang> |
<lang parigp>matqr(M)</lang> |
||
=={{header|Python}}== |
|||
{{libheader|numpy}} |
|||
<lang python>#!/usr/bin/env python3 |
|||
import numpy as np |
|||
def qr(A): |
|||
m, n = A.shape |
|||
Q = np.eye(m) |
|||
for i in range(n - (m == n)): |
|||
H = np.eye(m) |
|||
H[i:, i:] = make_householder(A[i:, i]) |
|||
Q = np.dot(Q, H) |
|||
A = np.dot(H, A) |
|||
return Q, A |
|||
def make_householder(a): |
|||
v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0])) |
|||
v[0] = 1 |
|||
H = np.eye(a.shape[0]) |
|||
H -= (2 / np.dot(v, v)) * np.dot(v[:, None], v[None, :]) |
|||
return H |
|||
# task 1: show qr decomp of wp example |
|||
a = np.array((( |
|||
(12, -51, 4), |
|||
( 6, 167, -68), |
|||
(-4, 24, -41), |
|||
)), dtype=np.float64) |
|||
q, r = qr(a) |
|||
print('q:\n', q.round(6)) |
|||
print('r:\n', r.round(6)) |
|||
# task 2: use qr decomp for polynomial regression example |
|||
def polyfit(x, y, n): |
|||
return lsqr(x[:, None]**np.arange(n + 1), y.T) |
|||
def lsqr(a, b): |
|||
q, r = qr(a) |
|||
_, n = r.shape |
|||
return np.linalg.solve(r[:n, :], np.dot(q.T, b)[:n]) |
|||
x = np.array((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), dtype=np.float64) |
|||
y = np.array((1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321), dtype=np.float64) |
|||
print('\npolyfit:\n', polyfit(x, y, 2))</lang> |
|||
{{out}} |
|||
<pre> |
|||
q: |
|||
[[-0.857143 0.394286 0.331429] |
|||
[-0.428571 -0.902857 -0.034286] |
|||
[ 0.285714 -0.171429 0.942857]] |
|||
r: |
|||
[[ -14. -21. 14.] |
|||
[ 0. -175. 70.] |
|||
[ 0. 0. -35.]] |
|||
polyfit: |
|||
[ 1. 2. 3.] |
|||
</pre> |
|||
=={{header|R}}== |
=={{header|R}}== |