diff --git a/linear_algebra/src/rayleigh_quotient_iteration.py b/linear_algebra/src/rayleigh_quotient_iteration.py
new file mode 100644
index 000000000000..081535efba01
--- /dev/null
+++ b/linear_algebra/src/rayleigh_quotient_iteration.py
@@ -0,0 +1,57 @@
+from typing import Tuple
+
+import numpy as np
+
+
+def rayleigh_quotient_iteration(
+    input_matrix: np.array,
+    vector: np.array,
+    value: float,
+    error_tol=1e-12,
+    max_iterations=100,
+) -> Tuple[float, np.array]:
+    """
+    https://en.wikipedia.org/wiki/Rayleigh_quotient_iteration
+
+    Rayleigh quotient iteration
+    Find an approximate eigenvector and eigenvalue of the input matrix. Which
+    eigenvalue, eigenvector pair is found depends on the vector and value
+    input. Has very rapid convergence.
+
+    >>> import numpy as np
+    >>> input_matrix = np.array([
+    ... [41,  4, 20],
+    ... [ 4, 26, 30],
+    ... [20, 30, 50]
+    ... ])
+    >>> vector = np.array([1, 1, 1])
+    >>> value = 30
+    >>> rayleigh_quotient_iteration(input_matrix, vector, value)
+    (33.60344292195729, array([-0.86306005,  0.45012643,  0.22915836]))
+    >>> value = 1
+    >>> rayleigh_quotient_iteration(input_matrix, vector, value)
+    (3.735693290153355, array([-0.23946876, -0.7641015 ,  0.59900218]))
+    """
+
+    N = len(input_matrix)
+
+    error = error_tol + 1
+    prev_eigenvalue = value
+    iterations = 0
+    while error > error_tol and iterations < max_iterations:
+        # Construct the iterating matrix
+        A = np.linalg.inv(input_matrix - np.eye(N) * value)
+        # Construct new eigenvector
+        vector = np.dot(A, vector)
+        # Normalize new eigenvector
+        vector = vector / np.linalg.norm(vector)
+
+        # Construct better approximation for eigenvalue
+        value = np.dot(np.conj(vector).T, np.dot(input_matrix, vector))
+
+        error = abs(value - prev_eigenvalue)
+        prev_eigenvalue = value
+
+        iterations += 1
+
+    return value, vector