Weighted Linear Least Squares - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

The Problem with Unweighted Residuals¶

As with most approximate methods, the moment we start to move away from an exact solution subtle effects start to show up!

Notice we are minimizing the residuals, but there is a subtle problem with the standard least squares formulation. Consider the inconsistent system from the pseudoinverse section:

\begin{aligned} 2 c + t &= 15 \\ c - t &= 10 \\ c + 2 t &= 14 \end{aligned}

(1)

If we were to multiply the second equation by 10, the mathematical line it represents remains geometrically identical. However, the residuals of the second equation would scale up by an order of magnitude:

\begin{aligned} 2 c + t &= 15 \\ 10 c - 10 t &= 100 \\ c + 2 t &= 14 \end{aligned}

(2)

Because the residuals for the second equation are now 10 times larger for the same geometric distance, the least squares algorithm will prioritize minimizing the error on the second equation, dragging our ‘center of mass’ solution towards it!

What does this remind you of?

Introducing Weights¶

Jacobi (diagonal) preconditioning!

P^{-1} = \begin{bmatrix} W_{1} & 0 & 0 \\ 0 & W_{2} & 0 \\ 0 & 0 & W_{3} \end{bmatrix}

(3)

If we define $r_1, r_2, r_3$ :

\begin{aligned} 2 c + t - 15 &= r_1 \\ c - t - 10 &= r_2 \\ c + 2 t - 14 &= r_3 \end{aligned}

(4)

We say the residuals can be weighted, i.e.: the least squares problem becomes,

\text{Min}_{\mathbf{x}} \ \sum W_i^2 r_i^2

(5)

Let’s code it and see how applying this diagonal matrix shifts the least squares solution!

A = np.array([[2, 1], [1, -1], [1, 2]])
b = np.array([15, 10, 14])

Pi = np.diag([1, 2, 1])
print("Preconditioner / Weights Matrix:\n", Pi)
A_weighted = Pi @ A
b_weighted = Pi @ b

x_lsq, _, _, _ = np.linalg.lstsq(A_weighted, b_weighted, rcond=None)
print("\nWeighted Least Squares Solution (c, t):", x_lsq)

x_unweighted, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
print("Unweighted Least Squares Solution (c, t):", x_unweighted)

# Prepare data for plotting
x_line = np.linspace(0, 10, 100)
y1 = 15 - 2*x_line
y2 = x_line - 10
y3 = (14 - x_line) / 2

# Plot the lines and the solutions
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_line, y=y1, mode='lines', name='2c + t = 15'))
fig.add_trace(go.Scatter(x=x_line, y=y2, mode='lines', name='c - t = 10'))
fig.add_trace(go.Scatter(x=x_line, y=y3, mode='lines', name='c + 2t = 14'))

fig.add_trace(go.Scatter(x=[x_unweighted[0]], y=[x_unweighted[1]], mode='markers', 
                         marker=dict(size=10, color='black'), 
                         name='Unweighted Solution'))

fig.add_trace(go.Scatter(x=[x_lsq[0]], y=[x_lsq[1]], mode='markers', 
                         marker=dict(size=14, color='red', symbol='star'), 
                         name='Weighted Solution (Equation 2 weighted)'))

fig.update_layout(title='Weighted vs Unweighted Least Squares for Inconsistent System', 
                  xaxis_title='c', yaxis_title='t',
                  yaxis=dict(range=[-10, 10]))
fig.show()

Preconditioner / Weights Matrix:
 [[1 0 0]
 [0 2 0]
 [0 0 1]]

Weighted Least Squares Solution (c, t): [9.33333333 0.33333333]
Unweighted Least Squares Solution (c, t): [8.33333333 1.33333333]

Pitfalls and opportunities¶

We have discovered that the magnitude of certain datapoints changes the best fit result, and awareness of this phenomenon is *critical for engineering applications.

Mitigation of these effects can be accomplished by sophisticated solver methods (beyond the scope of numpy / scipy), careful user guidance, or retroactive examination. However, the best practice are the ones you are already familiar with and centre around the removal of scale from data:

Using units
Scale transformations (linear vs log-scale)
Nondimensionalization

As usual, your engineering intuition and judgement will trump fancy numerics!

Measurement Uncertainty as Weights¶

A situtation where you will want to exploit weights is when you have redundant measurements and/or a known uncertainty.

If you have multiple measurements, $y_i^j$ , at a single point, $x_i$ , then you can find the mean,

\tilde{y}_i = \text{mean}_j (y_i^j)

(6)

and its variance,

\sigma_i^2 = \text{mean}\left([y_i^j - \tilde{y}_i]^2 \right)

(7)

and weight by the inverse of the variance, $W_i = \frac{1}{\sigma_i^2}$ .

Most package implementations of weighted least squares regression have an option to input the weights, but exact usage varies. Consult the documentation to determine how to use it in practice.