Exact rational least squares using Gram-Schmidt orthogonalization

qrls.py

from fractions import Fraction

def least_squares(A, b):
    b = list(map(Fraction, b))
    m = len(A)
    n = len(A[0]) if m else 0

    R = [[Fraction(0) for _ in range(n)] for _ in range(n)]

    Q_cols = []
    y = []
    pivots = []

    for j in range(n):
        v = [Fraction(A[i][j]) for i in range(m)]

        for i, q in enumerate(Q_cols):
            num = sum(qi * vi for qi, vi in zip(q, v))
            den = sum(qi * qi for qi in q)
            rij = num / den
            R[i][j] = rij
            v = [vi - rij * qi for qi, vi in zip(q, v)]

        if any(v_r != 0 for v_r in v):
            idx = len(Q_cols)
            pivots.append(j)
            Q_cols.append(v)

            num = sum(vi * bi for vi, bi in zip(v, b))
            den = sum(vi * vi for vi in v)
            y.append(num / den)

    rank = len(y)
    R = R[:rank]

    x = [Fraction(0) for _ in range(n)]
    for i in range(rank - 1, -1, -1):
        pivot = pivots[i]
        s = sum(R[i][k] * x[k] for k in range(pivot + 1, n))
        x[pivot] = y[i] - s

    return x