Andrew Adie

silhouette of trees under starry night sky

Least Squares

—

by

Andrew Adie
in

What is Least Squares?

Least squares is a mathematical (specifically statistical) method to identify a curve through data as an average of error through datapoints. It serves as a general basis for linear optimization problems.

History

Least squares began to surface as a result of a mathematical competition to describe the orbital patterns between Saturn and Jupiter. Leonard Euler, one of the most famous mathematicians in engineering mathematics, won the competition prize. His analysis included 8 unknown orbital elements, and 75 astrnomical observations. Although he intuitively understood the underlying discovery for least squares during his lifetime, the formal theory would not be fully developed until later.

The formal method was discovered in 1805.

Motivation

Suppose we have a system defined with the following relations:

x1+x2=4x_1 + x_2 = 4
2x1−x2=22x_1 – x_2 = 2
x1−2x2=1x_1 – 2x_2 = 1

This system has a definite analytical solution which would require 3 observations (4, 2, 1) on the right side of hte equalities. However, suppose our measurements are noisy and contain some error. How would these be quantified? We can separate the error from the “true” or “hidden” reading in the measurement

x1+x2=4+e1x_1 + x_2 = 4 +e_1
2x1−x2=2+e22x_1 -x_2=2 +e_2
x1−2x2=1+e3x_1 -2x_2 = 1 + e_3

where e’s represent error in the system. Now rearranging this system we have the following:

e1=x1+x2−4e_1 = x_1 + x_2 -4
e2=2x1−x2−2e_2 = 2x_1 -x_2-2
e3=x1−2x2−1e_3 = x_1 -2x_2 -1

Now the objective becomes to find the x variables so we minimize the errors:

minx1,x2⁡||e1+e2+e3||22\min_{x_1,x_2}||e_1+e_2+e_3||_2^2

This is the basis for the least squares method.

This method was justified by Gauss who used a maximum likeliehood distribution to do so.

Formulation of Least Squares Problem

suppose we have n unknown variables (X), and m observations (Y)

X=[x1,x2,...,xn]TX= [x_1, x_2, …, x_n]^T
Y=[y1,y2,...,ym]TY = [y_1, y_2,…,y_m]^T

We assume that the variables in X are related linearly. If this is the case, then we can more generally write:

Y=HX+WY = HX+W

Which encodes the relationships (H) between the true value (X) and the noise in our observations (W).

Note that these are all column vectors with the same dimension corresponding to the state variables.

Knowing we can’t directly measure the true value we need to estimate the variables. We use the following notation to indicate an estimator:

X^=[x^1,x^2,...,x^3]T\hat{X} = [\hat{x}_1, \hat{x}_2,…,\hat{x}_3]^T

since we defined error as the difference between the true and the estimated readings we can build the follwing relations:

e=Y−HX^e=Y-H\hat{X}

this gives us an idea that the error is the difference between the measurment and the estimator structure. In more advanced methods we assume that this is represented by (e=W), but for this derivation, we no longer need to make this assumption since we’re manipulating the error (e,W). So for each variable

e12+e22+...+em2=||e||22=||Y−HX^||22e_1^2 +e_2^2+…+e_m^2 = ||e||_2^2 = ||Y-H\hat{X}||_2^2

Returning to the opimization problem above, we can rewrite it as:

minX^⁡||e||22=minX^⁡||Y−HX^||22=minX^⁡J(X^)\min_{\hat{X}} ||e||_2^2 = \min_{\hat{X}}||Y-H\hat{X}||_2^2 = \min_{\hat{X}}J(\hat{X})

Where J denotes the objective function of the optimization problem.

Solution to Least Squares

Since this relationship between estimation variables is linear, we know the system is convex. As a result, if we can find a derivative of J equal to zero, we know we have found a global minimum of the system.

∇J(X^)=0\grad J(\hat{X}) = 0
∇J(X^)=∇||Y−HX^||22=2HTH−2X^Y=0\grad J(\hat{X}) = \grad ||Y-H\hat{X}||^2_2 = 2H^TH – 2\hat{X}Y=0

Therefore the optimality condition for least squares is:

HTH−HTY=0H^TH-H^TY=0

This result is a final solution if the observation matrix, H, is invertible. This means H must have at least n linearly independent rows.

Today, least squares may be done over an infinite-dimensional space, which will always have more unknowns than observations.

What happens when m < n?

Regularized Formulation

if m < n, then we regularize the problem to include a hyperparameter penalty. Our new formulation becomes:

minX^⁡||Y−HX^||22+λ||x^||22\min_{\hat{X}} ||Y-H\hat{X}||_2^2 +\lambda ||\hat{x}||_2^2

here we assume that lambda is a positive penalty parameter.

In this case the first order optimality condition becomes:

X^λ=(HTH+λI)−1HTY\hat{X}_{\lambda} = (H^TH+\lambda I)^{-1}H^TY

This is ALWAYS invertible from lambda times the identity when lambda is greater than zero.

The maginitude of necessary lambda is proportional to the inverse of the number of observations.

Evaluating the limit as lambda decreases to zero:

limλ→0⁡X^λ=X^0=H†Y\lim_{\lambda \to 0} \hat{X}_{\lambda} = \hat{X}_0 = H^{\dag}Y

where H-dagger is the psuedo-inverse of H.

This becomes teh solution to the following problem:

minX^⁡||X^||22s.t.||Y−HX^||22≤||Y−HX~||22∀X~\min_{\hat{X}} ||\hat{X}||_2^2 \, \, \, s.t. \, \, \, ||Y-H\hat{X}||_2^2 \leq ||Y-H \tilde{X}||_2^2 \, \, \, \forall \, \, \, \tilde{X}

meaning it is the globally minimum norm to this problem.

if lambda is a positive definite matrix, then you end up with the R matrix from LQR instead.

(λ||X||22→XTRX)(\lambda ||X||_2^2 \to X^TRX)

Comments

Leave a Reply

Discover more from Andrew Adie

Subscribe now to keep reading and get access to the full archive.

Continue reading