TheAlgorithms-Python

1
import numpy as np
2

3

4
def qr_householder(a: np.ndarray):
5
    """Return a QR-decomposition of the matrix A using Householder reflection.
6

7
    The QR-decomposition decomposes the matrix A of shape (m, n) into an
8
    orthogonal matrix Q of shape (m, m) and an upper triangular matrix R of
9
    shape (m, n).  Note that the matrix A does not have to be square.  This
10
    method of decomposing A uses the Householder reflection, which is
11
    numerically stable and of complexity O(n^3).
12

13
    https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections
14

15
    Arguments:
16
    A -- a numpy.ndarray of shape (m, n)
17

18
    Note: several optimizations can be made for numeric efficiency, but this is
19
    intended to demonstrate how it would be represented in a mathematics
20
    textbook.  In cases where efficiency is particularly important, an optimized
21
    version from BLAS should be used.
22

23
    >>> A = np.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)
24
    >>> Q, R = qr_householder(A)
25

26
    >>> # check that the decomposition is correct
27
    >>> np.allclose(Q@R, A)
28
    True
29

30
    >>> # check that Q is orthogonal
31
    >>> np.allclose(Q@Q.T, np.eye(A.shape[0]))
32
    True
33
    >>> np.allclose(Q.T@Q, np.eye(A.shape[0]))
34
    True
35

36
    >>> # check that R is upper triangular
37
    >>> np.allclose(np.triu(R), R)
38
    True
39
    """
40
    m, n = a.shape
41
    t = min(m, n)
42
    q = np.eye(m)
43
    r = a.copy()
44

45
    for k in range(t - 1):
46
        # select a column of modified matrix A':
47
        x = r[k:, [k]]
48
        # construct first basis vector
49
        e1 = np.zeros_like(x)
50
        e1[0] = 1.0
51
        # determine scaling factor
52
        alpha = np.linalg.norm(x)
53
        # construct vector v for Householder reflection
54
        v = x + np.sign(x[0]) * alpha * e1
55
        v /= np.linalg.norm(v)
56

57
        # construct the Householder matrix
58
        q_k = np.eye(m - k) - 2.0 * v @ v.T
59
        # pad with ones and zeros as necessary
60
        q_k = np.block([[np.eye(k), np.zeros((k, m - k))], [np.zeros((m - k, k)), q_k]])
61

62
        q = q @ q_k.T
63
        r = q_k @ r
64

65
    return q, r
66

67

68
if __name__ == "__main__":
69
    import doctest
70

71
    doctest.testmod()
72