gmuraru · breezko · Oct 4, 2020 · gmuraru · Oct 6, 2020 · gmuraru
diff --git a/...hon/nice_implementations/test/__init__.py → ...s/Python/nice_implementations/__init__.py b/...hon/nice_implementations/test/__init__.py → ...s/Python/nice_implementations/__init__.py
diff --git a/Languages/Python/nice_implementations/env.yml b/Languages/Python/nice_implementations/env.yml
@@ -11,3 +11,5 @@ dependencies:
     - mypy==0.782
     - mypy-extensions==0.4.3
     - pylint==2.6.0
+    - numpy==1.19.2
+
diff --git a/...ice_implementations/test/trie/__init__.py → ...plementations/src/householder/__init__.py b/...ice_implementations/test/trie/__init__.py → ...plementations/src/householder/__init__.py
diff --git a/Languages/Python/nice_implementations/src/householder/householder.py b/Languages/Python/nice_implementations/src/householder/householder.py
@@ -0,0 +1,80 @@
+
+
+class Householder:
+    def __init__(self, x, y):
+        """
+        Initialise a householder transformation which transforms the vector x to
+        point in the direction of a vector y, by reflecting x in a hyperplane
+        through the origin, whose normal is denoted by v. The matrix of this
+        linear transformation is denoted by P and is equal to I -
+        2*v@v.T/(v.T@v) == I - beta*v@v.T, where beta = 2/(v.T@v). By
+        definition, P@x == alpha*y for some scalar alpha.
+
+        According to the derivation in section 5.1.2 of Matrix Computations by
+        Golub and Van Loan (4th edition), expanding P = I - 2*v@v.T/(v.T@v) and
+        P@x == alpha*y implies that v is parallel to x - alpha*y, and the
+        magnitude of v is not unique, so we can set v = x - alpha*y. Expanding
+        P@x - alpha*y == 0 implies that alpha**2 == (x.T@x)/(y.T@y)
+        """
+        alpha = np.linalg.norm(x) / np.linalg.norm(y)
+        self.v = x - alpha * y
+        denom = self.v.T @ self.v
+        beta = 0 if denom == 0 else 2.0 / denom
+        self.beta_times_v = beta * self.v
+
+    def get_matrix(self):
+        """
+        Return the Householder matrix P == I - beta*(v @ v.T)
+        """
+        return np.identity(self.v.size) - np.outer(self.beta_times_v, self.v)
+
+    def pre_mult(self, A):
+        """
+        Compute the matrix product P @ A:
+
+        P @ A == (I - beta * v @ v.T) @ A == A - (beta * v) @ (v.T @ A)
+        """
+        return A - np.outer(self.beta_times_v, self.v.T @ A).reshape(A.shape)
+
+    def __matmul__(self, A):
+        return self.pre_mult(A)
+
+
+def qr(A, verbose=False):
+    """
+    Compute the QR decomposition, A = Q @ R, where Q @ Q.T = I, and R is upper
+    triangular, using Householder transformations.
+    """
+    m, n = A.shape
+    y = np.zeros(m)
+    y[0] = 1
+    R = A.copy()
+    k = min(m, n)
+    h_list = [None] * k
+    # Iterate through columns of A
+    for i in range(k):
+        # Apply Householder transform to sub-matrix to remove sub-diagonal
+        if verbose:
+            print("i={}\n".format(i), R, end="\n\n")
+        h = Householder(R[i:, i], y[: m - i])
+        h_list[i] = h
+        R[i:, i:] = h @ R[i:, i:]
+
+    # Create Q matrix from list of Householder transformations
+    Q = np.identity(A.shape[0])
+    for i in reversed(range(k)):
+        Q[i:, i:] = h_list[i] @ Q[i:, i:]
+
+    if verbose:
+        print(
+            "A=", A, "R=", R, "Q=", Q, "Errors=", qr_errors(A, Q, R), sep="\n\n"
+        )
+
+    return Q, R
+
+
+def qr_errors(A, Q, R):
+    recon_error = np.max(np.abs(A - Q @ R))
+    ortho_error = np.max(np.abs(np.identity(A.shape[0]) - Q @ Q.T))
+    triang_error = np.max(np.abs([e for i in range(R.shape[1]) for e in R[i + 1 :, i]]))
+    return recon_error, ortho_error, triang_error
diff --git a/Languages/Python/nice_implementations/tests/Small QR example.txt b/Languages/Python/nice_implementations/tests/Small QR example.txt
@@ -0,0 +1,17 @@
+A=
+
+[[ 12 -51   4]
+ [  6 167 -68]
+ [ -4  24 -41]]
+
+Q=
+
+[[ 0.8571428571 -0.3942857143 -0.3314285714]
+ [ 0.4285714286  0.9028571429  0.0342857143]
+ [-0.2857142857  0.1714285714 -0.9428571429]]
+
+R=
+
+[[ 14  21 -14]
+ [  0 175 -70]
+ [  0   0  35]]
diff --git a/Languages/Python/nice_implementations/tests/__init__.py b/Languages/Python/nice_implementations/tests/__init__.py
diff --git a/Languages/Python/nice_implementations/tests/householder/__init__.py b/Languages/Python/nice_implementations/tests/householder/__init__.py
diff --git a/Languages/Python/nice_implementations/tests/householder/test_householder.py b/Languages/Python/nice_implementations/tests/householder/test_householder.py
@@ -0,0 +1,100 @@
+""" Module to test the Householder class for Householder transformations """
+
+from ...src.householder.householder import Householder
+import numpy as np
+from numpy.linalg import norm
+
+np.random.seed(0)
+TOL = 1e-15
+N = 10
+M = 12
+
+
+def gaussian_vector(m):
+    """ Generate an m-dimensional Gaussian vector """
+    return np.random.normal(size=m)
+
+
+def gaussian_matrix(m, n, r=None):
+    """
+    Generate an m*n matrix using Gaussian random numbers. If r is None, the
+    matrix will almost surely be full rank; otherwise, if r is a non-negative
+    integer, the matrix will have rank r.
+    """
+    if r is None:
+        return np.random.normal(size=[m, n])
+    else:
+        A = np.zeros([m, n])
+        for _ in range(r):
+            A += np.outer(gaussian_vector(m), gaussian_vector(n))
+
+        return A
+
+
+def test_symmetric():
+    """ Test that the Householder matrix is symmetric """
+    x, y = gaussian_vector(N), gaussian_vector(N)
+    h = Householder(x, y)
+    P = h.get_matrix()
+    assert np.allclose(P, P.T, 0, TOL)
+
+
+def test_orthogonal():
+    """ Test that the Householder matrix is orthogonal """
+    x, y = gaussian_vector(N), gaussian_vector(N)
+    h = Householder(x, y)
+    P = h.get_matrix()
+    assert np.allclose(P @ P.T, np.identity(N), 0, TOL)
+
+
+def test_self_inverse():
+    """ Test that the Householder transformation is self-inverse """
+    x, y = gaussian_vector(N), gaussian_vector(N)
+    h = Householder(x, y)
+    assert np.allclose(x, h @ (h @ x), 0, TOL)
+    P = h.get_matrix()
+    assert np.allclose(P @ P, np.identity(N), 0, TOL)
+
+
+def test_householder_multiplication():
+    """
+    Test that the efficient Householder transform implementation is equivalent
+    to matrix multiplication
+    """
+    x, y = gaussian_vector(N), gaussian_vector(N)
+    h = Householder(x, y)
+    z = gaussian_matrix(N, M)
+    assert np.allclose(h @ z, h.get_matrix() @ z, 0, TOL)
+
+
+def test_householder_direction():
+    """
+    Test that the Householder transform points vectors in the correct direction
+    """
+    x, y = gaussian_vector(N), gaussian_vector(N)
+    h = Householder(x, y)
+    z = h @ x
+    assert np.allclose(z / norm(z), y / norm(y), 0, TOL)
+
+
+def test_zero_input():
+    """
+    Verify that if x is a zero-vector, then the Householder transform for any
+    given y is the identity transformation
+    """
+    x, y = np.zeros(N), gaussian_vector(N)
+    h = Householder(x, y)
+    assert np.allclose(h.get_matrix(), np.identity(N), 0, TOL)
+
+
+def test_zero_out_column():
+    """
+    Test using a Householder transform to zero out the sub-diagonal of the first
+    column of a matrix
+    """
+    A = gaussian_matrix(M, N)
+    y = np.zeros(M)
+    y[0] = 1
+    h = Householder(A[:, 0], y)
+    B = h @ A
+    assert np.allclose(B[1:, 0], 0, 0, TOL)
diff --git a/Languages/Python/nice_implementations/tests/householder/test_qr.py b/Languages/Python/nice_implementations/tests/householder/test_qr.py
@@ -0,0 +1,69 @@
+""" Module to test the QR decomposition using Householder transformations """
+
+from ...src.householder.householder import qr, qr_errors
+from .test_householder import gaussian_matrix
+import numpy as np
+
+np.random.seed(0)
+np.set_printoptions(precision=10, linewidth=1000, suppress=True, threshold=1000)
+TOL = 1e-13
+N = 13
+M = 17
+
+
+def _test_qr(m, n):
+    """
+    Test that the QR decomposition works for m*n matrices, including both
+    full-rank and low-rank matrics
+    """
+    A = gaussian_matrix(m, n)
+    Q, R = qr(A)
+    for e in qr_errors(A, Q, R):
+        assert e < TOL
+
+    for rank in range(4):
+        A = gaussian_matrix(m, n, rank)
+        Q, R = qr(A)
+        for e in qr_errors(A, Q, R):
+            assert e < TOL
+
+
+def test_qr_square():
+    """
+    Test that the QR decomposition works for square matrices, including both
+    full-rank and low-rank matrics
+    """
+    _test_qr(M, M)
+
+
+def test_qr_skinny():
+    """
+    Test that the QR decomposition works for skinny matrices (more rows than
+    columns), including both full-rank and low-rank matrics
+    """
+    _test_qr(M, N)
+
+
+def test_qr_fat():
+    """
+    Test that the QR decomposition works for fat matrices (more columns than
+    rows), including both full-rank and low-rank matrics
+    """
+    _test_qr(N, M)
+
+
+def test_small_example():
+    """ Test small 3x3 example from the Wikipedia page on QR decomposition """
+    A = np.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]])
+    Q, R = qr(A)
+    print("A=", A, "Q=", Q, "R=", R, sep="\n\n", file=open("Small QR example.txt", "w"))
+    for e in qr_errors(A, Q, R):
+        assert e < TOL
+
+
+def test_integer_matrix():
+    """ Test randomly generated integer matrix, and print output to file """
+    A = np.random.choice(4, [10, 14]).astype(np.float)
+    Q, R = qr(A, verbose=True)
+    for e in qr_errors(A, Q, R):
+        assert e < TOL