Why do we want a matrix?

Writing ten equations separately is inconvenient. Instead:

Xβ ≈ y

This is the central transformation of Milestone 2.

Why do we add a column of 1s?

Our equation is Sales = β₀ + β₁·P + β₂·G. We want to write it as y = Xβ. So we construct:

        X = 
[ 1   P₁   G₁ ]
[ 1   P₂   G₂ ]
[ ...  ...  ... ]
[ 1   P₁₀  G₁₀ ]

The first row [1, P₁, G₁] multiplied by [β₀, β₁, β₂]ᵀ gives β₀ + β₁P₁ + β₂G₁. The column of ones isn't a programming trick — it's how we mathematically encode the intercept.

Your actual design matrix

Your X has shape (10, 3):

• 10 rows = one per observation/year
• 3 columns: intercept, Population, GDP

What is β?

We collect all unknown coefficients into one vector:

β = [β₀, β₁, β₂]ᵀ

This is called the parameter vector. The entire regression problem is essentially: Find β.

What is y?

Your target vector contains the 10 Sales values: [4.63, 4.72, 4.83, 5.04, 5.25, 5.40, 5.56, 6.51, 7.04, 7.22]ᵀ. It has shape (10,).

Why do we write ≈ instead of =?

Because our model won't perfectly predict every observation. For example, ŷ₂₀₂₀ = 5.50 while y₂₀₂₀ = 5.56. Therefore Xβ ≠ y in general. Instead, Xβ is as close as possible to y under our least-squares criterion.

The fundamental problem

You have 10 equations and only 3 unknowns. This is an overdetermined system. We generally cannot satisfy all ten equations exactly.

Why can't we simply solve Xβ = y?

A lot of introductory linear algebra teaches Ax = b → x = A⁻¹b. But that assumes A is a square matrix. Your X is (10×3) — it isn't square. You cannot simply do X⁻¹.

How do we turn it into something solvable?

We start with Xβ ≈ y and the least-squares objective. The mathematical derivation eventually gives:

XᵀXβ = Xᵀy

These are called the normal equations.

Why does XᵀX become square?

Your X is (10×3). Transpose Xᵀ is (3×10). Multiply: (3×10)(10×3) = (3×3). That's square! And Xᵀy is (3×10)(10×1) = (3×1).

So the normal equations have the structure: (3×3)(3×1) = (3×1). Now we have a standard square linear system.

But why multiply by Xᵀ at all?

This isn't done merely because "we need a square matrix." It naturally emerges from minimizing the least-squares error. XᵀX isn't an arbitrary trick — it's the mathematical consequence of taking derivatives and setting them to zero.

You'll see the actual derivation in Milestone 6. For now, remember:

Least Squares → Normal Equations → XᵀXβ = Xᵀy

The complete transformation

CSV DATA (10 years)
        ↓
DESIGN MATRIX X  (10 × 3)
        ↓
XᵀX and Xᵀy
        ↓
A = XᵀX  (3 × 3)
b = Xᵀy  (3 × 1)
        ↓
Aβ = b   →  MILESTONE 3

What is a transpose?

Transpose swaps rows ↔ columns. If X is (10×3), then Xᵀ is (3×10).

What does @ mean?

In Python, X @ beta means matrix multiplication. So X.T @ X means XᵀX.

What does hstack do?

np.hstack([A, b]) places A and b side-by-side to create the augmented matrix [A|b]. For a 3×3 A and 3×1 b, the result is 3×4.

Scale problem

Your features operate on wildly different scales: Population ≈ 320,000,000 vs GDP ≈ 60,000. When you calculate XᵀX, Population gets squared to ~10¹⁷. This creates a numerical stability problem — your condition number eventually shows ~1.9×10²¹. That's the problem Milestone 6 will fix.

Elementary row operations

We can perform three operations without changing the solution:

1. Swap rows
2. Multiply a row by a nonzero number
3. Add a multiple of one row to another

The basic idea

Gauss-Jordan tries to transform the augmented matrix [A|b] into [I|β], where I is the identity matrix. Then the answer is simply the final column.

Partial pivoting

For each column, we find the largest absolute value and swap that row to the pivot position. This protects against numerical instability. If the pivot is essentially zero, the algorithm raises a singularity error.

Normalize and eliminate

Divide the pivot row by the pivot value (making it 1), then eliminate all other entries in that column (both above and below — that's what makes it Gauss-Jordan, not just Gaussian elimination).

Why is pivoting important?

Real computers use floating-point numbers. Dividing by a tiny number can amplify numerical errors. Choosing a larger pivot generally makes elimination more stable.

The idea

Instead of directly solving Aβ = b, decompose:

A = LU

where L is lower triangular and U is upper triangular.

Two simpler problems

1. Solve Lz = Pb using forward substitution
2. Solve Uβ = z using backward substitution

P is a permutation matrix representing row swaps from pivoting. We must apply the same swaps to b: pb = P @ b.

Why is LU useful for repeated solves?

If A doesn't change but b does, you decompose A=LU once and reuse it for every new b. That's much more efficient than re-doing the entire elimination.

Doolittle's method

Constructs L with ones on the diagonal. The elimination factors are stored in L while U gets the resulting upper-triangular matrix. Think of L as a record of the elimination operations.

Forward and backward substitution

Forward substitution solves Lz = Pb starting from the first equation downward. Backward substitution solves Uβ = z starting from the last equation upward.

The library reference

np.linalg.solve(A, b) acts as your reference implementation. If all three solvers agree:

β_GJ ≈ β_LU ≈ β_NumPy

That's strong evidence your implementation is correct.

Singularity and near-singularity

A singular matrix doesn't have a unique solution — the equations are redundant. A near-singular matrix has equations that are almost redundant, making tiny input errors cause large changes in the answer.

The condition number

Your raw XᵀX shows a condition number of ~1.9×10²¹. That's enormous. After standardization in Milestone 6, it drops to ~22. That's one of the most important results in your project.

The basic idea

Give me an interval containing the root, and I'll keep cutting that interval in half.

The bracket condition

If f(a) and f(b) have opposite signs (f(a)·f(b) < 0), there must be a root between them — this relies on the Intermediate Value Theorem.

The algorithm

1. Compute midpoint m = (a+b)/2
2. Evaluate f(m)
3. If f(a)·f(m) < 0, root is in [a, m]; otherwise in [m, b]
4. Repeat until the interval is small enough

Why is Bisection "safe"?

It never leaves the original bracket. Every subsequent interval remains inside it. Bisection can't suddenly produce an absurd answer. That's why we describe it as:

Slow but robust

Convergence rate

The interval is halved every iteration. After n iterations, the interval width is (b−a)/2ⁿ. This is linear convergence — the number of correct digits increases at a steady but slow rate.

A radically different idea

Instead of shrinking a safe interval, Newton says: Look at the slope of the function and jump toward the root.

The Newton formula

x_new = x − f(x)/f'(x)

This is one of the most important formulas in numerical methods.

Where does it come from?

Use the tangent-line approximation. Set y = 0 and solve for x_new. The derivative f'(x) gives the slope, and the formula tells you where the tangent line crosses zero.

For our problem

Because our model is linear in GDP, the derivative is simply β₂ (the GDP coefficient). Each iteration asks: how far are we from the target, and how strongly does Sales respond to GDP?

Why did Newton take only 3 iterations?

For a perfectly linear function, the tangent line is the function itself. Therefore Newton can jump directly to the root.

Bisection = safety · Newton = speed

Newton's weakness

Newton is not guaranteed to stay near the root. For complicated functions, it can converge to a different root, oscillate, or diverge. It only sees the function locally — it doesn't know the global shape.

Why both methods gave the same answer

Your project obtained ≈$87,138.33 from both. They agree because your model is linear and well-behaved. If the function were nonlinear with multiple roots, the methods could behave very differently.

The deeper lesson

Milestone 3 gave us the model parameters β. Milestone 4 turns those parameters into a decision question:

What input is required to achieve a desired output?

The clever idea

Build a special polynomial for each data point that gives that point a weight of 1 and all other points a weight of 0. Then combine them:

p(x) = y₀L₀(x) + y₁L₁(x) + y₂L₂(x) + ...

At x = x₀: p(x₀) = y₀(1) + y₁(0) + y₂(0) = y₀. The polynomial passes through every point.

The basis formula

Lᵢ(x) = ∏ⱼ≠ᵢ (x − xⱼ)/(xᵢ − xⱼ)

A tiny example

Points: (1,2) and (3,6). At x=2: L₀(2) = (2−3)/(1−3) = 0.5, L₁(2) = (2−1)/(3−1) = 0.5. So p(2) = 2(0.5) + 6(0.5) = 4.

The idea

Instead of constructing all basis polynomials directly, calculate divided differences — generalized slopes, then changes in slopes, then changes in those changes.

First divided difference = slope

f[x₀, x₁] = (y₁ − y₀)/(x₁ − x₀)

Second divided difference

f[x₀, x₁, x₂] = (f[x₁,x₂] − f[x₀,x₁])/(x₂ − x₀)

Conceptually: how is the slope changing? This corresponds to curvature.

Newton's formula

p(x) = f[x₀] + f[x₀,x₁](x−x₀) + f[x₀,x₁,x₂](x−x₀)(x−x₁) + ...

The key advantage: you can extend the polynomial by adding another term when a new data point arrives.

Same polynomial

Given n distinct x-values, there is exactly one polynomial of degree at most n−1 that passes through all n points. Lagrange and Newton are just different representations of the same polynomial.

The dangerous part

You have 10 observations. Using all 10 produces a degree-9 polynomial. It passes through every point perfectly — but between and beyond the points, it can behave wildly.

This is called Runge's phenomenon.

Your project demonstrates this

Full degree-9 polynomial at GDP=66,000 gives ≈$3.9677T. Local (k=4) interpolation gives ≈$4.9780T. That's a huge difference — and 66,000 is inside the data range.

Extrapolation disaster

At GDP=86,032 (just outside the range), the degree-9 polynomial produces ≈$24,478.83T. That's absurd — but the computer correctly evaluated the polynomial. The problem is the model.

A mathematically valid model can still be practically useless.

The solution: local interpolation

Instead of using all 10 points, use the k nearest observations. With k=4, you get a degree-3 polynomial — much more restrained. The idea: nearby information is usually more relevant.

The central lesson

Exact fit ≠ good general behavior

This is one of the most important ML principles hiding inside numerical analysis.

The immediate problem

We have 10 equations but only 3 unknowns (β₀, β₁, β₂). We usually cannot make all 10 equations exactly true.

What does "error" mean?

For observation i: the residual is eᵢ = yᵢ − ŷᵢ. This is how wrong the model was for that observation.

Why not simply add all errors?

Errors can cancel: +2 + (−2) = 0 looks like no error. So we square them:

SSE = Σ(yᵢ − ŷᵢ)²

Our optimization problem: find β₀, β₁, β₂ that minimize SSE.

The connection to linear algebra

The design matrix X and target vector y give us:

ŷ = Xβ

The error is y − Xβ. Least squares minimizes |y − Xβ|².

Setting the gradient to zero

We take partial derivatives with respect to each β and set them to zero. This gives the normal equations:

XᵀXβ = Xᵀy

These are called the normal equations — one of the most important equations in classical linear regression.

The big connection

Minimize squared error
        ↓ (calculus)
XᵀXβ = Xᵀy
        ↓ (linear algebra)
Aβ = b
        ↓ (Milestone 3 solvers)
Best β
        ↓
Predictions

The scale problem

Your two features have wildly different scales: Population ≈ 3×10⁸ vs GDP ≈ 6×10⁴. That's roughly a 10,000× difference. When you calculate XᵀX, the scale difference gets squared to ~10⁸ or more.

Standardization

Convert each feature to a z-score: z = (x − μ) / σ. After standardization, both features are roughly centered around 0 with standard deviation 1.

The dramatic improvement

Raw condition number: 1.8999 × 10²¹ (extremely ill-conditioned)
Standardized condition number: 22.0227 (dramatically better)

This is one of the most important results in the project.

The fix was data preprocessing, not a better solver

When we encounter numerical problems, the better response is often: "Is the system poorly conditioned because of how we represented the data?" Standardization changes the geometry and numerical conditioning of the optimization problem.

Prediction uses the same scaler

During prediction, we reuse the same mean and standard deviation learned during training. We do not recalculate them. This principle is extremely important in ML: Fit preprocessing on training data; reuse it on new data.

How well does the model explain the data?

R² = 1 − SS_res / SS_tot

where SS_res is unexplained squared error and SS_tot is total variation around the mean.

What does R² mean?

Your result: R² = 0.9755. This means approximately 97.55% of the variation in Sales is explained by the fitted linear relationship with Population and GDP.

But be careful: don't say "Population and GDP cause 97.55% of Sales." R² measures explained variation by the model, not causation.

Don't overinterpret with n=10

With only 10 observations, R² = 0.9755 does not prove the model will predict future Sales with 97.55% accuracy. It only describes fit to these observations.

Standardized coefficients

Because the model standardizes features, β₁ and β₂ correspond to: change in predicted Sales associated with a one-standard-deviation increase in that feature, holding the other constant. For raw-space sensitivity, use the chain rule: ∂Sales/∂x = β/σ.

Milestone 5: exact interpolation vs Milestone 6: optimized approximation

Partial derivatives as sensitivities

For our linear model:

∂Sales/∂Population = β₁
∂Sales/∂GDP = β₂

Think of β₁ and β₂ as sensitivities — how strongly the output reacts to each input.

The chain rule for standardization

Because the model was trained on standardized features, we need to convert back to raw units:

∂Sales/∂Population = β₁/σ_Population
∂Sales/∂GDP = β₂/σ_GDP

The key equation

σ_Sales ≈ √((∂Sales/∂P · σ_P)² + (∂Sales/∂G · σ_G)²)

This is the delta method / first-order error propagation formula. We square independent uncertainties and combine them using root-sum-of-squares.

The idea

Instead of doing calculus, generate thousands of plausible Population/GDP values and see how much the prediction changes.

Simulation 1: Pop=334.8M, GDP=$85,431 → Sales=7.48T
Simulation 2: Pop=335.2M, GDP=$84,721 → Sales=7.42T
Simulation 3: Pop=334.9M, GDP=$85,177 → Sales=7.45T
...
Simulation 20,000: → Sales=7.44T

The standard deviation of all those predictions gives us our empirical estimate of Sales uncertainty.

Why Monte Carlo?

The analytical method assumes linearity and independence. Monte Carlo doesn't depend on those assumptions — it tests the uncertainty empirically. If both methods agree, that gives us confidence.

Your results

Analytical uncertainty ≈ 0.08333T · Monte Carlo uncertainty ≈ 0.08448T. They're extremely close — good validation.

Variance contributions

Population contribution ≈ 0.000012 · GDP contribution ≈ 0.006932

GDP contributes roughly 580× more variance to predicted Sales than Population.

What is an uncertainty budget?

Like a financial budget, but for uncertainty. You want to know: where is most of the uncertainty coming from?

Total prediction uncertainty
├── Population → tiny contribution
└── GDP → dominant contribution

If your objective is to improve the reliability of Sales prediction: improving the quality of the GDP input is potentially much more valuable than improving Population.

The volume knob analogy

Population is a small knob — move it a little, Sales barely changes. GDP is a large knob — move it the same relative amount, Sales changes much more. If the GDP knob itself is shaky, your final prediction becomes shaky too.

We use Year as the input

Your data shows Population generally increasing and GDP generally increasing over time. So we can fit: Population = a + b(Year) and GDP = c + d(Year).

Why use regression again?

This is a beautiful connection. Milestone 6 already taught us how to fit y = β₀ + β₁x. Now we reuse the same machinery. We're not inventing a new algorithm — we're reusing what we already built.

Why use degree 1 (linear trend)?

You only have 10 years of data. A high-degree polynomial could memorize historical fluctuations rather than learn the underlying trend. A simple linear trend is a reasonable baseline. Remember Milestone 5: more flexibility isn't automatically better.

The most important terminology

Interpolation — predict inside known data

known ─── unknown ─── known
          ↑
       interpolation

Generally safer.

Extrapolation — predict outside known data

known ───────────────→ unknown
                       ↑
                  extrapolation

Generally riskier.

Milestone 5 vs Milestone 8

Milestone 5 used interpolation (estimating inside the observed GDP range). Milestone 8 uses extrapolation (predicting beyond 2023). That's why we should be much more cautious.

Why extrapolation is dangerous

Inside the range, neighboring observations constrain the curve. Outside the range, there are no observations telling the curve what to do. The historical data cannot tell you with certainty what happens after the last observation.

Forecasting is not just extending a line. It is extending an assumption.

Why we care about residuals

Historical Population trend errors (residuals) give us a rough measure of how uncertain our trend prediction is. We calculate their standard deviation: pop_resid_std and gdp_resid_std.

Uncertainty grows with forecast distance

The code uses: σ_future = σ_historical × (1 + 0.5d), where d = years beyond the last observation. For 2025 (d=2): uncertainty is 2× historical. For 2026 (d=3): 2.5×.

This is a naive uncertainty-growth heuristic — a modeling choice, not a fundamental law. Sophisticated models (ARIMA, Bayesian forecasting, etc.) could estimate uncertainty differently.

The complete pipeline

MILESTONE 1: Define the relationship
        ↓
MILESTONE 2: Matrix construction
        ↓
MILESTONE 3: Solve equations
        ↓
MILESTONE 6: Fit Sales model
        ↓
Sales = f(Pop, GDP)
        ↓
   ┌────┴────┐
   ↓         ↓
M4: Target  M7: Uncertainty
GDP?        propagation
   ↓         ↓
   └────┬────┘
        ↓
MILESTONE 8
        ↓
Future Year → Pop forecast → GDP forecast
        ↓
   Sales model → Future Sales
        ↓
   Uncertainty interval

The single most important equation

Year → {Population, GDP} → Sales

Milestone 6 was: Population, GDP → Sales. Milestone 8 adds: Year → Population, GDP. We had the Sales model already. The missing problem was: where do future Population and GDP come from?