Table of Contents

Topic 8. Linear Regression and Least Squares¶

Author: Brady Anderson
    b.anderson@byu.edu

Introduction¶

Linear Regression, or the Least Squares problem, manifests itself quite often in engineering work. First, an example:

A physical (real-world) engineering process needs to be modeled. Now, having the physical system, it is possible to apply known inputs to the process "plant," and measure the corresponding outputs. (Assuming some type of output measurements are possible with some sensor). Linear Regression allows plant parameters to be deduced given a collected set of inputs and outputs. These plant parameters become your model, whether for simulation or for an estimator.

The idea is to find linear coefficients for a set of linearly independent vectors such that the combination of these coefficients and vectors (the estimated data) comes as close as possible to a true dataset.

Linear Regression takes advantage of properties of the Projection Theorem and the induced norms on $l_{2}$ and $L_{2}$ in order to compute the aforementioned parameters (or whatever the solution may be to your use case). The key fact to note is that the error of the estimated data is orthogonal to the true data.

Explanation of the theory¶

Give a detailed discussion (i.e., equations galore) of the topic.  The emphasis here is clarity for future students learning the topic.

First a review of the Projection Theorem (Thm 2.9 in Moon):

Let $S$ be a Hilbert space and let $V$ be a closed subspace of $S$. For any vector $\mathbf{x} \in S$, $\exists$ a unique vector $\mathbf{v}_0 \in V$ closest to $\mathbf{x}$; that is, $\|\mathbf{x}-\mathbf{v_0}\| \le\|\mathbf{x}-\mathbf{v}\| \ \forall \ \mathbf{v} \ \in V$. Furthermore, the point $\mathbf{v}_0$ is the minimizer of $\|\mathbf{x}-\mathbf{v}\|$ if and only if $\mathbf{x}-\mathbf{v_0} \in V^{\perp}$.

If you are unfamiliar or uncomfortable with this definition, please see the section on projection operators. Also, section 2.14 of the book covers the theorem well.

Say you have a normed, linear vector space ($S,\|\cdot\|$), and a matrix of linearly independent vectors $T=[\mathbf{p}_1,\mathbf{p}_2,\dots,\mathbf{p}_m]$ with $V=span(T) \in S$. In linear regression, we seek the coeffients $c_i$ such that

$\mathbf{\hat{x}}=T\mathbf{\hat{c}}$

according to the constraint that $\mathbf{\hat{c}}=\underset{c}{\operatorname{argmin}}\|\mathbf{x}-T\mathbf{c}\|$.

This will minimize the total error between the true vector $\mathbf{x}$ and the estimate $\mathbf{\hat{x}}$:
$\|\mathbf{e}\|=\|\mathbf{x}-\mathbf{\hat{x}}\|$.
It is geometrically observable that the error is minimized when the error is orthogonal to $V$ when the $l_{2}$ and $L_{2}$ norms are used.

Now, if $\mathbf{x} \in V$, it is possible to have $\|\mathbf{e}\|=0$. However, in general, (especially in engineering), this is not the case.

When we use the induced norm for $\|\mathbf{e}\|$, we can express the minimization in terms of the orthogonality condition, using the Projection Theorem: the minimum-norm error must be orthogonal to each vector $\mathbf{p}_j$

$\langle \mathbf{x}-\sum_{i=1}^{m}{c_i\mathbf{p}_i}, \mathbf{p}_j \rangle = 0, \ \ \ \ \ j=1,2,\dots,m$

We can write these in what are known as the normal equations, with $m$ equations in $m$ unknowns:

$$\begin{bmatrix} \langle \mathbf{p}_1, \mathbf{p}_1 \rangle & \cdots & \langle \mathbf{p}_m, \mathbf{p}_1 \rangle \\ \vdots & \ddots & \vdots \\ \langle \mathbf{p}_1, \mathbf{p}_m \rangle & \cdots & \langle \mathbf{p}_m, \mathbf{p}_m \rangle \end{bmatrix} \begin{bmatrix} c_1 \\ \vdots \\ c_m \end{bmatrix} = \begin{bmatrix} \langle \mathbf{x}, \mathbf{p}_1 \rangle \\ \vdots \\ \langle \mathbf{x}, \mathbf{p}_m \rangle \end{bmatrix}$$

The left matrix is known as the Grammian of the set T, the set of vectors composing the T matrix. It usually is denoted by the letter $R$: $R\mathbf{c}=\mathbf{p}_x$

where $\mathbf{p}_x$ is the cross-correlation vector.

Due to the properties of the inner product, and that $R_{ij}=\langle \mathbf{p}_j,\mathbf{p}_i \rangle$, the Grammian matrix is Hermitian symmetric: $R^{H}=R$. (The Hermitian of a matrix means its conjugate transpose.)

To solve the normal equations, $R$ must be invertible. Recall that positive-definite matrices are always invertible. This leads to $\mathbf{Thm 3.1}$:

A Grammian matrix $R$ is always positive-semidefinite (that is, $\mathbf{x}^{H}R\mathbf{x} \ge 0 \ \forall \ \mathbf{x} \in \mathbb{C}^{m}$). It is positive-definite if and only if the vectors $\mathbf{p}_1,\dots,\mathbf{p}_m$ are linearly independent.

Proof:

Let $\mathbf{y}=[y_1,\dots,y_m]^{T}$ be an arbitrary vector. Then

$$\mathbf{y}^{H}R\mathbf{y}=\sum_{i=1}^{m}{\sum_{j=1}^{m}{\bar{y_i}y_j \langle \mathbf{p}_i, \mathbf{p}_j \rangle}} = \langle \sum_{j=1}^{m}{y_j\mathbf{p}_j}, \sum_{i=1}^{m}{y_i\mathbf{p}_i} \rangle = \left\| \sum_{j=1}^{m}{y_j\mathbf{p}_j} \right\|^{2} \ge 0$$

Hence $R$ is positive-semidefinite.

If $R$ is not positive-semidefinite, then there is a nonzero vector $\mathbf{y}$ such that

$\mathbf{y}^{H}R\mathbf{y}=0$

so that

$\sum_{i=1}^{m}{y_i\mathbf{p}_i}=0$;

thus, the $\mathbf{p}_i$ are linearly independent.

Conversely, if $R$ is positive-definite, then

$\mathbf{y}^{H}R\mathbf{y} \gt 0$

for all nonzero $\mathbf{y}$ and by the first line equation of this proof,

$\sum_{i=1}^{m}{y_i\mathbf{p}_i} \not= 0$

This means that the $\mathbf{p}_i$ are linearly independent.     QED

As a result of this theorem, the Grammian is invertible if all the vectors $\mathbf{p}_i$ are linearly independent.

As an extension, if the set of vectors $\mathbf{p}_i$ are orthogonal, then the Grammian is diagonal, significantly reducing the amount of computation required to find the coefficients. They are simply obtained by:

$c_{j}=\frac{\langle \mathbf{x}, \mathbf{p}_j \rangle}{\langle \mathbf{p}_j, \mathbf{p}_j \rangle}$

The orthogonality principle for least-squares is now formalized with Thm 3.2:

Let $\mathbf{p}_1,\dots,\mathbf{p}_m$ be data vectors in a vector space $S$. Let $\mathbf{x}$ be any vector in $S$. In the representation

$\mathbf{x} = \sum_{c=1}^{m}{c_{i}p_{i}} + \mathbf{e} = \mathbf{\hat{x}} + \mathbf{e}$

the induced norm of the error $\|\mathbf{e}\|$ is minimized when the error $\mathbf{e}=\mathbf{x}-\mathbf{\hat{x}}$ is orthogonal to each of the data vectors,

$\langle \mathbf{x}-\sum_{i=1}^{m}{c_{i}\mathbf{p}_{i}}, \mathbf{p}_j \rangle=0, \ \ \ \ \ j=1,2,\dots,m$

Proof (via Cauchy-Schwarz inequality):

In the case that $\mathbf{x} \in span(\mathbf{p}_1, \dots, \mathbf{p}_m)$, the error is zero and hence is orthogonal to the data vectors. This case is therefor trivial and is excluded from what follows.

If $\mathbf{x} \notin span(\mathbf{p}_1, \dots, \mathbf{p}_m)$, let $\mathbf{y}$ be a fixed vector that is orthogonal to all of the data vectors,

$\langle \mathbf{y},\mathbf{p}_i \rangle=0 \ \ \ \ \ i=1,\dots,m$

such that

$\mathbf{x}=\sum_{i=1}^{m}{a_i\mathbf{p}_i}+\mathbf{y}$

for some set of coefficients ${a_1, \dots, a_m}$. Let $\mathbf{e}$ be a vector satisfying

$\mathbf{x}=\sum_{i=1}^{m}{c_i\mathbf{p}_i}+\mathbf{e} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ eq(1)$

for some set of coefficients ${c_1, \dots, c_m}$. Then by the Cauchy-Schwarz inequality,

$$\begin{eqnarray} \|\mathbf{e}\|^2\|\mathbf{y}\|^2 &\ge& | \langle \mathbf{e},\mathbf{y} \rangle |^2 \\ &=&\left| \langle \mathbf{x},\mathbf{y} \rangle - \langle \sum_{c=1}^{m}{c_{i}p_{i}}, \mathbf{y} \rangle \right|^2 \\ &=&| \langle \mathbf{x},\mathbf{y} \rangle |^2 \end{eqnarray}$$

The lower bound is independent of the coefficients ${c_i}$, and hence no set of coefficients can make the bound smaller. By the quality condition for the Cauchy-Schwarz inequality, the lower bound is achieved --implying the minimum $\|\mathbf{e}\|$ -- when

$\mathbf{e}=\alpha\mathbf{y}$

for some scalar $\alpha$. Since $\mathbf{e}$ must satisfy eq(1), it must be the case that:

$\alpha=1$

$a_i=c_i$

$\mathbf{e}=\mathbf{y}$

hence the error is orthogonal to the data.     QED

It is important to note that because $\mathbf{\hat{x}}$ is a linear combination of the data vectors, it is orthogonal to the error vector:

$\langle \mathbf{\hat{x}}, \mathbf{e} \rangle = 0$

Or, in other words, the residual vector (errors from each point) is orthogonal to the column space. See demonstration below.

Simple Numerical Examples¶

Provide some simple python code and examples that emphasize the basic concepts.

Note Figure 3.8(b), provided in the book.

It may be confusing that the minimized error is the vertical distance betweent the regressed line and the data points. The following example will show how these vertical distances does not contradict the earlier statements of orthogonal error minimization.

Let data point vector $\mathbf{y}=[1,4,5]^T$, measured at time $\mathbf{t}=[0,1,2]^T$ be approximated by a linear regression:

$\mathbf{y}=a_1\mathbf{t}+a_0$

This may be expressed as:

$$\begin{bmatrix} 1\\4\\5 \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & 1 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} a_1 \\ a_0 \end{bmatrix}$$

Here we have $\mathbf{p}_1=\mathbf{t}$ and $\mathbf{p}_2=[1,1,1]^T$

To compute the Grammian matrix, R, we need the inner product of $\mathbf{p}_1$ and $\mathbf{p}_2$ to themselves, and to each other.

We also compute the cross-correlation vector $\mathbf{p}_x$:

$$\mathbf{p}_x=[\langle \mathbf{y}, \mathbf{p}_1 \rangle, \langle \mathbf{y}, \mathbf{p}_2 \rangle]^T$$

Finally, we verify the linear independence of $\mathbf{p}_1, \mathbf{p}_2$ (thus guaranteeing inertibility of R) by asserting

$\langle \mathbf{\hat{p}}_1, \mathbf{\hat{p}}_2 \rangle \not= 1$

where $\mathbf{\hat{p}}_i$ is the unit vector in the direction of $\mathbf{p}_i$

Our p-vectors are in $\mathbb{R}^3$. When we say that our error is minimized by being orthogonal to the data vectors we mean that the residual is orthogonal to the column space of the matrix created by the data. In our case, this is the matrix composed of $\mathbf{p}_1$ and $\mathbf{p}_2$. When this is drawn on the 2D plot of the data, it looks like vertical distances from the true data to the regressed line. But if we plot the data vectors in 3D, in our case a plane, we can see that the residual error is perpendicular to this plane, shown by

$\langle \mathbf{p}_1, \mathbf{e} \rangle = 0$

$\langle \mathbf{p}_2, \mathbf{e} \rangle = 0$

$\mathbf{y}$ and $\mathbf{t}$, as well as the regressed line, are drawn below.

In the 3d plot, you can see the plane spanned by ${\mathbf{p}_1,\mathbf{p}_2}$, and that the estimated $\mathbf{\hat{y}}$ lies in this plane. You can also see the error vector is orthogonal to this plane created by the data vectors. Finally, if you add the error vector to the estimated vector, you get back the original, true data

$\mathbf{\hat{y}} + \mathbf{e} = \mathbf{y}$

seen by the fact the yellow vector representing the true data is out of the plane.

You can grab and rotate the plot with the mouse, to see these facts.

Finally, we make one more observation, specific to the matrix case.

Recall that the Grammian matrix, $R$, is the left hand matrix in the normal equations:

$$\begin{bmatrix} \langle \mathbf{p}_1, \mathbf{p}_1 \rangle & \cdots & \langle \mathbf{p}_m, \mathbf{p}_1 \rangle \\ \vdots & \ddots & \vdots \\ \langle \mathbf{p}_1, \mathbf{p}_m \rangle & \cdots & \langle \mathbf{p}_m, \mathbf{p}_m \rangle \end{bmatrix} \begin{bmatrix} c_1 \\ \vdots \\ c_m \end{bmatrix} = \begin{bmatrix} \langle \mathbf{x}, \mathbf{p}_1 \rangle \\ \vdots \\ \langle \mathbf{x}, \mathbf{p}_m \rangle \end{bmatrix}$$

If we denote our data matrix as $A$, note that $R$ may also be represented as $A^HA$. Additionally, we may present the right hand vector, $\mathbf{p}_x$, as $A^H\mathbf{x}$.

Also recall that if the columns in the data matrix, $A$, are linearly independent, then $R$ is invertible. Hence,

$R\mathbf{c}=\mathbf{p}_x$

$\Rightarrow A^HA\mathbf{c}=A^H\mathbf{x}$

$\Rightarrow \mathbf{c} = (A^HA)^{-1}A^H\mathbf{x}$

Below we show that the Moore-Penrose pseudo-inverse, $(A^HA)^{-1}A^H$, is equivalent to the process above when your data can be expressed as a matrix (i.e. not in function space).

An Engineering Application¶

Provide a more sophisticated example showing one engineering example of the topic, complete with python code.

Given a mass, spring, damper system,

$m\ddot{x} - b\dot{x} - kx = u$

with known inputs and outputs, estimate the system parameters, $m, k, b$.

Note: Outputs are usually measured with instruments that introduce noise into the measurements. So we call this an estimate because if we had perfect measurements, you could retrieve the exact values. Due to the noise in the measurements, and imperfect actuator, noise will be added into the data you use to retrieve the parameters.

It is easiest to take advantage of the matrix form of the Moore-Penrose Pseudo-Inverse:

$A\mathbf{c} = \mathbf{x}$

$\mathbf{c} = (A^HA)^{-1}A^H\mathbf{x}$

where, for our problem,

$$A=\begin{bmatrix} \mathbf{\ddot{x}} & \mathbf{-\dot{x}} & -\mathbf{x} \end{bmatrix}$$ $$c=\begin{bmatrix} \mathbf{m} \\ \mathbf{b} \\ \mathbf{k} \end{bmatrix}$$ $$\mathbf{x}=\mathbf{u}$$

Suggestion: Play with the sigma (standard deviation) multiplier $sigma_{mult}$ for the noise, as well as $t_{fin}$, to see limits of this method, and how longer data sampling times can help with the estimation.

Application to Estimating Simulation Results¶

Author: Jonathan Barrett
    jab128@byu.edu

Suppose you have created a model that evaluates the probability that an ADS-B message will be decoded successfully given a density of UAS/km$^2$. You have collected data for UAS densities from nearly 0 UAS/km$^2$ up to 5 UAS/km$^2$. Not only is the data noisy from randomization, but now the simulations are too resource intensive to continue running at higher densities.

You realize that you need to perform a least squares linear regression on the data to fit a line to the data, and to estimate the probability of a successful decode at 10 UAS/km$^2$. Use the following csv file, and the python template below: https://drive.google.com/file/d/1vBTGxYauh1itxIlzHmcZsaQXcVefohVv/view?usp=sharing

If your estimated decode probability seems wrong, then you may want to use a weighted least squares solution that follows the later trends in the data most closely.

Congratulations! You found coefficients to form a linear equation that best fits the data, but you can actually fit the data better by finding coefficients for polynomials of different orders. Re-write your program so that it can fit the data to any order polynomial that you choose.

Once you have done so, set the polynomial order to 10 and run the program. Does the curve seem to fit the data better? Why is that?

Check the approximations for the ADS-B successful decode probability at 10 UAS/km$^2$. Do they seem better or worse than the approximations from the linear equation? Why is that?

Now set the polynomial order to 100 and run the program. Does the curve seem to fit the data well?