Python/maths/qr_decomposition.py
Caeden Perelli-Harris 490e645ed3
Fix minor typing errors in maths/ (#8959)
* updating DIRECTORY.md

* types(maths): Fix pylance issues in maths

* reset(vsc): Reset settings changes

* Update maths/jaccard_similarity.py

Co-authored-by: Tianyi Zheng <tianyizheng02@gmail.com>

* revert(erosion_operation): Revert erosion_operation

* test(jaccard_similarity): Add doctest to test alternative_union

* types(newton_raphson): Add typehints to func bodies

---------

Co-authored-by: github-actions <${GITHUB_ACTOR}@users.noreply.github.com>
Co-authored-by: Tianyi Zheng <tianyizheng02@gmail.com>
2023-08-15 14:27:41 -07:00

72 lines
2.1 KiB
Python

import numpy as np
def qr_householder(a: np.ndarray):
"""Return a QR-decomposition of the matrix A using Householder reflection.
The QR-decomposition decomposes the matrix A of shape (m, n) into an
orthogonal matrix Q of shape (m, m) and an upper triangular matrix R of
shape (m, n). Note that the matrix A does not have to be square. This
method of decomposing A uses the Householder reflection, which is
numerically stable and of complexity O(n^3).
https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections
Arguments:
A -- a numpy.ndarray of shape (m, n)
Note: several optimizations can be made for numeric efficiency, but this is
intended to demonstrate how it would be represented in a mathematics
textbook. In cases where efficiency is particularly important, an optimized
version from BLAS should be used.
>>> A = np.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]], dtype=float)
>>> Q, R = qr_householder(A)
>>> # check that the decomposition is correct
>>> np.allclose(Q@R, A)
True
>>> # check that Q is orthogonal
>>> np.allclose(Q@Q.T, np.eye(A.shape[0]))
True
>>> np.allclose(Q.T@Q, np.eye(A.shape[0]))
True
>>> # check that R is upper triangular
>>> np.allclose(np.triu(R), R)
True
"""
m, n = a.shape
t = min(m, n)
q = np.eye(m)
r = a.copy()
for k in range(t - 1):
# select a column of modified matrix A':
x = r[k:, [k]]
# construct first basis vector
e1 = np.zeros_like(x)
e1[0] = 1.0
# determine scaling factor
alpha = np.linalg.norm(x)
# construct vector v for Householder reflection
v = x + np.sign(x[0]) * alpha * e1
v /= np.linalg.norm(v)
# construct the Householder matrix
q_k = np.eye(m - k) - 2.0 * v @ v.T
# pad with ones and zeros as necessary
q_k = np.block([[np.eye(k), np.zeros((k, m - k))], [np.zeros((m - k, k)), q_k]])
q = q @ q_k.T
r = q_k @ r
return q, r
if __name__ == "__main__":
import doctest
doctest.testmod()