mirror of
https://github.com/TheAlgorithms/Python.git
synced 2024-10-06 13:49:30 +00:00
07e991d553
* ci(pre-commit): Add pep8-naming to `pre-commit` hooks (#7038) * refactor: Fix naming conventions (#7038) * Update arithmetic_analysis/lu_decomposition.py Co-authored-by: Christian Clauss <cclauss@me.com> * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * refactor(lu_decomposition): Replace `NDArray` with `ArrayLike` (#7038) * chore: Fix naming conventions in doctests (#7038) * fix: Temporarily disable project euler problem 104 (#7069) * chore: Fix naming conventions in doctests (#7038) Co-authored-by: Christian Clauss <cclauss@me.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
72 lines
2.1 KiB
Python
72 lines
2.1 KiB
Python
import numpy as np
|
|
|
|
|
|
def qr_householder(a):
|
|
"""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()
|