graduate_code/solutions/numerics.py

"""G05: educational Householder least squares for bounded full-rank matrices."""
import math


def least_squares(A, b, rtol=1e-12):
    m, n = len(A), len(A[0]) if A else 0
    if not 1 <= n <= 64 or not n <= m <= 10000 or len(b) != m:
        raise ValueError("require m >= n, 1..64 columns, <=10000 rows")
    if any(len(row) != n for row in A):
        raise ValueError("ragged matrix")
    if not math.isfinite(rtol) or not 0 < rtol < 1:
        raise ValueError("relative rank threshold must lie in (0,1)")
    if any(not math.isfinite(x) for row in A for x in row) or any(not math.isfinite(x) for x in b):
        raise ValueError("finite matrix and target required")
    scale = max(abs(x) for row in A for x in row)
    if scale == 0:
        raise ValueError("rank-deficient zero matrix")
    R = [[x/scale for x in row] for row in A]
    target = [x/scale for x in b]
    if any(not math.isfinite(x) for x in target):
        raise ValueError("target exceeds supported dynamic range")
    for j in range(n):
        z = [R[i][j] for i in range(j,m)]
        norm = math.hypot(*z)
        if norm <= rtol:
            raise ValueError("rank deficiency at stated relative threshold")
        alpha = -math.copysign(norm,z[0])
        v = list(z)
        v[0] -= alpha
        vnorm = math.hypot(*v)
        v = [x/vnorm for x in v]
        for k in range(j,n):
            projection = 2*math.fsum(v[i-j]*R[i][k] for i in range(j,m))
            for i in range(j,m):
                R[i][k] -= projection*v[i-j]
        projection = 2*math.fsum(v[i-j]*target[i] for i in range(j,m))
        for i in range(j,m):
            target[i] -= projection*v[i-j]
        R[j][j] = alpha
        for i in range(j+1,m):
            R[i][j] = 0.
    x = [0.]*n
    for j in range(n-1,-1,-1):
        remainder = math.fsum(R[j][k]*x[k] for k in range(j+1,n))
        x[j] = (target[j]-remainder)/R[j][j]
    if not all(math.isfinite(value) for value in x):
        raise ValueError("nonfinite coefficients")
    return x


if __name__ == "__main__":
    print(least_squares([[1,0],[1,1],[1,2]],[1,2,2]))