Turn Several Measurements Into Several Predictions

Botie has two measurements for every greenhouse. A simple illustrative score gives each measurement a coefficient, adds its contributions, then adds an offset. We can calculate one greenhouse's score by hand before doing the same for a whole table.

This is where the distinction between elementwise multiplication and matrix multiplication matters. We will keep the observations, feature columns, coefficients, and output shapes visible at every step.

Temperature contribution = temperature × 2
Humidity contribution = humidity × 0.1
Score = temperature contribution + humidity contribution + 3

Choose one greenhouse row and predict what will happen if one coefficient increases. The calculation below combines the contributions for every row using the same scoring rule.

Elementwise multiplication combines matching entries without adding them. A dot product multiplies matching entries and then adds those products. For one observation, the feature and coefficient vectors must have equal lengths.

The coefficients in this example are supplied solely to illustrate the calculation. They are not a trained rule or a physical greenhouse control law.

import numpy as np

features = np.array([10.,40.])
weights = np.array([2.,0.1])
contributions = features * weights
print(contributions.tolist())
print(np.dot(features, weights))
print(np.dot(features, weights) + 3)

[20.0, 4.0]
24.0
27.0

The @ operator performs matrix multiplication. For a two-dimensional feature table and a one-dimensional coefficient vector, it computes one dot product per row. The feature count must match the coefficient count.

For matrix times matrix, the left matrix's column count must match the right matrix's row count. The surviving outer dimensions describe the output.

import numpy as np

X = np.array([[10.,40.],[20.,60.]])
weights = np.array([2.,0.1])
scores = X @ weights + 3
print(scores.tolist())
print(scores.shape)
print((X @ np.eye(2)).tolist())

[27.0, 49.0]
(2,)
[[10.0, 40.0], [20.0, 60.0]]

A norm measures the size of a vector. The Euclidean norm squares its coordinates, adds them, and takes the square root. np.linalg.norm performs this calculation; naming an axis lets us calculate it separately for each row.

When coordinates are different kinds of measurement, a raw distance can be dominated by their units. Later scaling lessons explain how to make such comparisons meaningful.

import numpy as np

errors = np.array([[3.,4.],[0.,5.]])
print(np.linalg.norm(errors, axis=1).tolist())
print(np.sqrt(np.square(errors).sum(axis=1)).tolist())

[5.0, 5.0]
[5.0, 5.0]

If a square coefficient matrix has an independent constraint for each unknown, solve can find values satisfying that system. Do not explicitly compute an inverse merely to apply it to one right-hand side. A singular matrix does not supply a unique solution.

With more observations than unknowns, exact agreement may be impossible. lstsq finds coefficients minimising the sum of squared residuals. It also returns information about the fit and the rank, the number of independent directions supported by the matrix.

import numpy as np

A = np.array([[1.,1.],[1.,-1.]])
b = np.array([30.,10.])
settings = np.linalg.solve(A, b)
print(settings.tolist())
print((A @ settings).tolist())

[20.0, 10.0]
[30.0, 10.0]

import numpy as np

x = np.array([0.,1.,2.])
y = np.array([1.,2.,2.])
X = np.column_stack([np.ones_like(x), x])
weights, residual_sums, rank, singular_values = np.linalg.lstsq(X, y, rcond=None)
print(np.round(weights, 4).tolist())
print(np.round(X @ weights, 4).tolist())
print(rank)

[1.1667, 0.5]
[1.1667, 1.6667, 2.1667]
2

An exact solve asks for coefficients that satisfy every equation. What if the measurements contain noise and the equations disagree? Then we can look for predictions whose squared residuals have the smallest total. That is the least-squares problem.

We arrange one observation per row of the design matrix, with a column for each coefficient. A column of ones allows the model to fit an intercept. Once the solver returns its coefficients, let's use them to calculate the predictions and residuals. Those show how the fitted rule behaves on our actual measurements.

import numpy as np
x = np.array([0., 1., 2.])
y = np.array([1., 2., 5.])
X = np.column_stack([np.ones(x.size), x])
w, _, rank, _ = np.linalg.lstsq(X, y, rcond=None)
print(np.round(w, 3).tolist())
print(np.round(y - X @ w, 3).tolist())
print(int(rank))

[0.667, 2.0]
[0.333, -0.667, 0.333]
2

If two feature columns are identical, changing one coefficient and compensating with the other can leave every prediction unchanged. A solver can still return a useful prediction rule, but the individual coefficients are not uniquely determined by those rows.

Inspect rank and singular values when this distinction matters. More decimal places do not recover missing information. A regularisation rule can choose among possible coefficients, but that is an additional modelling choice. It does not make the data identify two independent effects.