An eigenvector of a square matrix A is a direction that A only stretches: A @ v = lam v, with the eigenvalue lam the stretch. Multiply any vector by A over and over, normalising each time, and it turns toward the eigenvector with the largest eigenvalue. This is power iteration.

import numpy as np

A = np.array([[2.0, 1.0], [1.0, 2.0]])
v = np.ones(2) + np.array([1.0, 0.0])
for _ in range(50):
    v = A @ v
    v = v / np.linalg.norm(v)
print(v.round(3), (v @ A @ v).round(3))
print(np.linalg.eigvalsh(A))  # [1. 3.]

Any matrix, square or not, factors as U @ diag(S) @ Vt, with orthonormal U and Vt and singular values S from largest to smallest. Keeping only the first k gives the best rank-k approximation, and the error left over is the root of the sum of the squared dropped singular values. Low-rank fine-tuning methods such as LoRA rely on weight updates being well approximated this way.

Write power_iteration(A, steps=200) for a symmetric matrix, starting from np.ones, returning (value, vector): the largest eigenvalue as a float (as v @ A @ v) and its unit eigenvector. Write low_rank(M, k) returning the best rank-k approximation of M, using np.linalg.svd(M, full_matrices=False).

The explanation opens here once every check passes.