1import json, math
  2import numpy as np
  3
  4
  5def householder_qr_factors(A):
  6    """Unblocked Householder QR, returning reflector factors and R.
  7    The factors (v,tau) are an implicit representation of Q; no Q is stored.
  8    """
  9    A = np.array(A, dtype=np.float64, copy=True)
 10    m, n = A.shape
 11    p = min(m, n)
 12    factors = []
 13    for j in range(p):
 14        x = A[j:, j]
 15        nx = np.linalg.norm(x)
 16        if nx == 0:
 17            factors.append((np.zeros_like(x), 0.0)); continue
 18        # Sign choice avoids cancellation.
 19        alpha = -math.copysign(nx, x[0])
 20        v = x.copy(); v[0] -= alpha
 21        vv = float(v @ v)
 22        tau = 0.0 if vv == 0 else 2.0 / vv
 23        A[j:, j:] -= np.outer(v, tau * (v @ A[j:, j:]))
 24        A[j, j] = alpha
 25        A[j+1:, j] = 0.0
 26        factors.append((v, tau))
 27    return factors, A
 28
 29
 30def apply_factors_left(factors, X, transpose=False):
 31    """Apply product of stored reflectors to X (reflectors are symmetric)."""
 32    Y = np.array(X, dtype=np.float64, copy=True)
 33    seq = factors if transpose else factors[::-1]
 34    for j, (v, tau) in enumerate(seq):
 35        row = j if transpose else len(factors)-1-j
 36        # In reverse traversal, the factor's original row is row.
 37        if transpose:
 38            row = j
 39        Y[row:] -= np.outer(v, tau * (v @ Y[row:]))
 40    return Y
 41
 42
 43def householder_block_basis(Y):
 44    fac, R = householder_qr_factors(Y)
 45    # Q(:,1:b) is obtained by applying Q to first b coordinate vectors.
 46    b = Y.shape[1]
 47    E = np.eye(Y.shape[0], b)
 48    Q = apply_factors_left(fac, E, transpose=False)
 49    return Q, fac, R
 50
 51
 52def gs_block_basis(Y):
 53    """One-pass classical Gram-Schmidt, intentionally no reorthogonalization."""
 54    m, b = Y.shape; Q = np.zeros((m, b)); r = 0
 55    for j in range(b):
 56        v = Y[:, j].copy()
 57        if r: v -= Q[:, :r] @ (Q[:, :r].T @ v)
 58        nv = np.linalg.norm(v)
 59        if nv > 1e-14:
 60            Q[:, r] = v / nv; r += 1
 61    return Q[:, :r]
 62
 63
 64def adaptive(A, sigma, block, method="householder", kmax=None, seed=0):
 65    rng = np.random.default_rng(seed)
 66    m, n = A.shape; kmax = min(m, n) if kmax is None else min(kmax, m)
 67    R = A.copy(); initial = float(np.sum(R*R)); E = initial
 68    blocks = []; cols = 0
 69    while E > sigma*sigma and cols < kmax:
 70        b = min(block, kmax-cols)
 71        # Fresh Gaussian range sketch of the current residual.
 72        Y = R @ rng.standard_normal((n, b))
 73        if method == "householder": Q, _, _ = householder_block_basis(Y)
 74        else: Q = gs_block_basis(Y)
 75        if Q.shape[1] == 0: break
 76        # Q is already orthogonal to previous blocks because it is formed from R.
 77        R -= Q @ (Q.T @ R)
 78        blocks.append(Q); cols += Q.shape[1]; E = float(np.sum(R*R))
 79    Qall = np.concatenate(blocks, axis=1) if blocks else np.zeros((m,0))
 80    return {"Q": Qall, "residual_energy": E, "initial_energy": initial,
 81            "rank": Qall.shape[1], "orth_error": float(np.linalg.norm(Qall.T@Qall-np.eye(Qall.shape[1]))) }
 82
 83
 84def matrix_with_spectrum(m, n, decay):
 85    rng = np.random.default_rng(1234 + int(decay*1000))
 86    U,_=np.linalg.qr(rng.standard_normal((m,m))); V,_=np.linalg.qr(rng.standard_normal((n,n)))
 87    s=np.exp(-decay*np.arange(min(m,n)))
 88    return U[:, :min(m,n)] @ np.diag(s) @ V[:min(m,n), :]
 89
 90
 91def run():
 92    # Prediction 1: exact Householder orthogonality stays near epsilon as blocks grow;
 93    # one-pass GS worsens strongly for ill-conditioned panels.
 94    rng=np.random.default_rng(7); m=96; b=4
 95    cond_rows=[]
 96    for cond in [1e2,1e6,1e10,1e14]:
 97        x=np.linspace(0,1,m); Y=np.column_stack([x**j + 1e-14*rng.standard_normal(m) for j in range(b)])
 98        # Make the conditioning control explicit through singular values.
 99        U,_=np.linalg.qr(rng.standard_normal((m,b))); W,_=np.linalg.qr(rng.standard_normal((b,b)))
100        Y=U@np.diag(np.geomspace(1,1/cond,b))@W.T
101        qh,_,_=householder_block_basis(Y); qg=gs_block_basis(Y)
102        cond_rows.append([cond,float(np.linalg.norm(qh.T@qh-np.eye(b))),float(np.linalg.norm(qg.T@qg-np.eye(qg.shape[1])))])
103    # Prediction 2: for geometric spectrum, rank grows monotonically as sigma tightens.
104    A=matrix_with_spectrum(72,48,0.16); normA=np.linalg.norm(A,'fro')
105    tol_rows=[]
106    for rel in [0.5,0.25,0.12,0.06,0.03]:
107        out=adaptive(A, rel*normA, 4, "householder", seed=11)
108        tol_rows.append([rel,out['rank'],math.sqrt(out['residual_energy'])/normA, out['orth_error']])
109    # Prediction 3: captured residual is below prescribed sigma (up to sampling/roundoff).
110    checks=[]
111    for rel in [0.4,0.2,0.1]:
112        out=adaptive(A, rel*normA, 4, "householder", seed=21)
113        checks.append([rel, math.sqrt(out['residual_energy'])/normA, out['rank']])
114    # Secondary baseline comparison at same max rank.
115    base=adaptive(A, .1*normA, 4, "gs", seed=21); idea=adaptive(A,.1*normA,4,"householder",seed=21)
116    result={"predicted": {"orthogonality": "Householder O(eps), GS grows with cond", "rank": "tighter sigma -> nondecreasing rank", "residual": "relative residual <= sigma/||A||"}, "orthogonality_sweep":cond_rows,"tolerance_sweep":tol_rows,"residual_checks":checks,"comparison":{"gs":{k:base[k] for k in ['rank','residual_energy','orth_error']},"householder":{k:idea[k] for k in ['rank','residual_energy','orth_error']}}}
117    print(json.dumps(result, indent=2))
118
119if __name__ == '__main__': run()