Inspect a poorly conditioned least-squares problem

Build a polynomial design matrix, inspect singular values, compare solvers and separate coefficient stability from prediction error.

Mathematics & statistics · Lesson 03🔋 Offline core analysis⏱ 50 minNumerical teaching model📦 numpy · pandas · matplotlib · IPython

Goal and how to run

This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.

Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.

Setup and parameters

In [1] · Python

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='ma03-conditioning'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

Load the embedded snapshot

In [2] · Python

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['models']

Read the method functions

Data and model boundary

Analytic functions and Monte Carlo points are mathematical teaching models, not experimental measurements. When Iris observations are used, row exchangeability and the sampling design remain assumptions; convenience records alone do not establish population representativeness. Numerical discretization error is distinct from sampling uncertainty.

Steps

1. Build a deterministic model benchmark

Use 25 x values on [0,1] and y=exp(x). A degree-10 Vandermonde design has powers with very different scales. These are analytic teaching values; no measurement noise or statistical confidence interval is implied.

In [3] · Python

inputs=snapshot("models")
DEGREE=10;x=np.linspace(0,1,25);y=np.exp(x)
design=np.vander(x,DEGREE+1,increasing=True)
singular=np.linalg.svd(design,compute_uv=False);condition=np.linalg.cond(design)
print("Design shape:",design.shape,"; condition number:",condition)
fig,ax=plt.subplots(figsize=(7,4));ax.semilogy(np.arange(1,len(singular)+1),singular,"o-")
ax.set(xlabel="Singular-value index",ylabel="Singular value",title="Different resolvability of polynomial directions");plt.tight_layout();plt.show()

Design shape: (25, 11) ; condition number: 21503150.614924222

2. Solve through SVD and compare normal equations

The SVD least-squares solver avoids squaring the condition number. The normal-equation comparison explicitly solves AᵀA; it can be unstable or singular even when A has full mathematical rank. Report solver failures instead of silently replacing them.

In [4] · Python

stable=np.linalg.lstsq(design,y,rcond=None)[0]
try:
 normal=np.linalg.solve(design.T@design,design.T@y);normal_failure=None
except np.linalg.LinAlgError as error:
 normal=None;normal_failure=str(error)
print("Condition(A):",condition,"; condition(A.T @ A):",np.linalg.cond(design.T@design))
print("Normal-equation failure:",normal_failure)
dense=np.linspace(0,1,501);dense_design=np.vander(dense,DEGREE+1,increasing=True)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(dense,np.exp(dense),color="black",label="Analytic exp(x)")
ax.plot(dense,dense_design@stable,"--",label="SVD least squares")
if normal is not None:ax.plot(dense,dense_design@normal,":",label="Normal equations")
ax.scatter(x,y,s=12);ax.set(xlabel="x",ylabel="y",title="Predictions can agree despite coefficient instability");ax.legend();plt.tight_layout();plt.show()

Condition(A): 21503150.614924222 ; condition(A.T @ A): 462238754711098.0
Normal-equation failure: None

3. Perturb the input response at a recorded scale

Add a deterministic 10⁻⁸-scale sinusoidal perturbation. Compare relative changes in response and coefficients, then inspect prediction error on the same interpolation interval. Large coefficient changes need not translate directly into the same prediction change.

In [5] · Python

perturbed=y+1e-8*np.sin(17*x)
altered=np.linalg.lstsq(design,perturbed,rcond=None)[0]
relative_input=np.linalg.norm(perturbed-y)/np.linalg.norm(y)
relative_coefficients=np.linalg.norm(altered-stable)/np.linalg.norm(stable)
print("Relative response perturbation:",relative_input,"; relative coefficient change:",relative_coefficients)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(dense,dense_design@stable-np.exp(dense),label="Original response")
ax.plot(dense,dense_design@altered-np.exp(dense),label="Perturbed response")
ax.set(xlabel="x in interpolation interval",ylabel="Prediction − analytic value",title="Observable prediction errors");ax.legend();plt.tight_layout();plt.show()

Relative response perturbation: 3.900268533006143e-09 ; relative coefficient change: 0.007641567140337687

4. Check orthogonality and export diagnostics

A least-squares residual should be approximately orthogonal to design columns, subject to floating-point precision. Record the matrix condition, solver choice, coefficient vectors and observation/prediction errors separately. Do not interpret these deterministic coefficients as causal effects.

In [6] · Python

residual=design@stable-y
assert np.linalg.norm(design.T@residual)<1e-8
assert np.isfinite(stable).all() and np.max(np.abs(dense_design@stable-np.exp(dense)))<1e-5
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"power":np.arange(DEGREE+1),"SVD_coefficient":stable,"perturbed_coefficient":altered,"normal_equation_coefficient":normal if normal is not None else np.full(DEGREE+1,np.nan)}).to_csv(OUTPUT_DIR/"coefficients.csv",index=False)
pd.DataFrame({"x":dense,"analytic":np.exp(dense),"SVD_prediction":dense_design@stable,"perturbed_prediction":dense_design@altered}).to_csv(OUTPUT_DIR/"interpolation.csv",index=False)
(OUTPUT_DIR/"diagnostics.json").write_text(json.dumps({"degree":DEGREE,"condition_design":condition,"condition_gram":float(np.linalg.cond(design.T@design)),"relative_input_perturbation":relative_input,"relative_coefficient_change":relative_coefficients,"normal_equation_failure":normal_failure},indent=2))
print("Least-squares checks passed; exports:",OUTPUT_DIR)

Least-squares checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ma03-conditioning

Exercises and reference explanation

Why can normal equations worsen conditioning?
Why can coefficients be unstable while interpolation remains accurate?

Reference explanation. Forming AᵀA approximately squares the condition number. Correlated basis directions can cancel in predictions; coefficient and prediction stability are different targets.

Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.

Optional AI coding exercises

Prompt 1

Using NumPy/Matplotlib, compare degrees 3, 6 and 10 while keeping the same analytic observations and solver diagnostics.

Prompt 2

Using the same packages, add a centered/scaled polynomial basis and compare condition numbers and interpolation errors.

Human acceptance. Record failures explicitly, check normal-equation residuals and keep synthetic/model labels.

Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.

Optional online extension

In [7] · Python

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://docs.scipy.org/doc/scipy/tutorial/integrate.html',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged

Sources, snapshots and next steps

IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256 a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.
Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256 4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.

Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.