Linear Algebra Algorithms in TheAlgorithms/Python

LU decomposition is preferred over Gaussian elimination when the same coefficient matrix is solved against multiple right-hand-side vectors, because the factorization is paid once. The implementation follows Doolittle's algorithm: L has ones on the diagonal and zeros above it, while U has the elimination multipliers on and above the diagonal. The outer loop advances i from 0 to n-1. For each i, the inner j loop below i fills L entries by subtracting the dot product of previously computed L and U slices from the original entry and dividing by the U diagonal. After setting the L diagonal to 1, a second j loop fills the corresponding U row. The zero-pivot check on line 101 catches cases where no decomposition exists, which occurs when a leading principal minor is zero. The doctest for a 3x3 matrix confirms the L and U factors.

Jacobi iteration solves a strictly diagonally dominant system by repeatedly updating each variable using only the other variables' previous values. The naive implementation has a Python inner loop over variables, but the vectorized version (lines 158-165) eliminates it entirely. Before the iteration loop, the code precomputes three structures: denominator holds the diagonal entries, no_diagonals is the coefficient matrix with zeros on the diagonal (extracted via a boolean mask), and ind is an index array that maps row positions to off-diagonal column indices. Inside each iteration, np.take(init_val, ind) gathers current variable values for broadcasting, and a single NumPy sum expression computes all weighted off-diagonal contributions at once. Dividing by the diagonal gives the next approximation. The function requires a strictly diagonally dominant matrix, verified by strictly_diagonally_dominant before the loop.

Matrix inversion has a narrow valid domain: the matrix must be square and have a non-zero determinant. NumPy's linalg.inv raises LinAlgError when those conditions fail, but that name is specific to NumPy's numeric layer. This file wraps the call to catch LinAlgError and raise ValueError with the message "Matrix is not invertible". Callers who do not import NumPy can handle the error by type alone, and the message is a domain-level description rather than a NumPy internal signal. The function also converts input and output between plain Python lists and NumPy arrays, keeping the interface free of NumPy types. This wrapper pattern is common in scientific Python where a library-level function must present a stable interface independent of the numeric backend. The doctest documents both the happy path and the singular matrix case.

The conjugate gradient method exploits the structure of symmetric positive definite (SPD) systems to converge in at most n iterations without storing a full factorization. The four update equations in the loop implement the core algorithm. First, alpha is the step size along the current search direction, computed as the squared residual norm divided by the directional curvature. Second, the solution guess advances by alpha times p0. Third, the residual shrinks by the same step. Fourth, beta is the ratio of consecutive squared residual norms, weighting the previous search direction when building the next one. The matrix-vector product w = A @ p0 is computed once per iteration and reused for both alpha and the residual update, halving the matrix multiplications. The outer validation asserts the matrix is SPD before the loop, since CG is only valid for SPD systems.

Power iteration finds the largest eigenvalue and its corresponding eigenvector without computing the full eigendecomposition. The algorithm works because multiplying a matrix by a vector repeatedly amplifies the component along the dominant eigenvector faster than any other component. Dividing by the vector norm after each multiplication prevents the values from growing without bound and keeps the iteration numerically stable. The Rayleigh quotient v^T A v provides the eigenvalue estimate at each step: because v is normalized, the Rayleigh quotient equals the eigenvalue when v is the exact eigenvector. Convergence is checked by comparing the change in the Rayleigh quotient between iterations relative to its current magnitude. The comment notes the Rayleigh quotient is faster than usual because the normalization step is already done. The doctest for a 3x3 matrix returns eigenvalue 79.66 and the matching eigenvector.

The rank of a matrix is the number of linearly independent rows, and Gaussian elimination reveals it by counting pivot positions that survive the reduction. This implementation starts with rank = min(rows, columns), the maximum possible value, and decrements whenever a pivot column cannot be filled. When the diagonal element is non-zero, the pivot is usable and the code eliminates below it in the standard way. When it is zero, the code searches below for a non-zero element to swap in. If none exists, the column is all-zero at this step, meaning one row is linearly dependent; rank decrements and the algorithm moves on. The doctest confirms rank 2 for a matrix with dependent rows, rank 4 for a full-rank 4x4, and rank 0 for an empty row list.