Linear Least Squares Methods with R: An Algebraic Approach - Part I

Algebraic Approach Principles

The least-squares method is one of the most well-known linear optimization methods because of its flexibility. Furthermore, it gives a reasonable approximation of a given function. Among the diverse applications that it can be used for is Regression in statistical learning, Direct Linear Transformation methods in projective geometry, and so on. We will demonstrate the principles of least squares methods and implement examples in R in this article.

Consider the above figure of a simple linear system where the $\mathbf{x}$$\mathbf{x}$ variable is an input variable, $A$$A$ is a measurements matrix, and $\mathbf{y}$$\mathbf{y}$ is the output variable. That is, the model for the equation is given by

$$\begin{array}{}\text{(1)}& \mathbf{y}=A\mathbf{x},\end{array}$$

which in a more explicit notation from $\text{(1)}$$\eqref{eq:ref1}$ can be expressed as

$$\begin{array}{}\text{(2)}& \left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_i\end{array}\right]=\left[\begin{array}{cccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,i}\\ a_{2,1}& a_{2,2}& \cdots & a_{2,i}\\ ⋮& ⋮& \cdots & ⋮\\ a_{j,1}& a_{j,2}& \cdots & a_{j,i}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_i\end{array}\right]\end{array}$$

We want to find and estimation of the variable $\mathbf{x}$$\mathbf{x}$ given the matrix measurement $A$$A$ and the output $\mathbf{y}$$\mathbf{y}$. To do that, we consider the equivalent estimated model input and output given by

$$\begin{array}{}\text{(3)}& \hat{\mathbf{y}}=A\hat{\mathbf{x}}.\end{array}$$

We then define a minimization criterion. In this case we use a quadratic minimization criterion, since we guarantee to locate global minima, that is,

$$\begin{array}{}\text{(4)}& \underset{x}{min}\parallel \mathbf{y}-\hat{\mathbf{y}}{\parallel }^2=\underset{x}{min}\parallel \mathbf{\text{y}}-A\hat{\mathbf{x}}{\parallel }^2.\end{array}$$

In order to estimate the output of $\hat{\mathbf{y}}$$\mathbf{\hat y}$ as closely as possible to $\mathbf{y}$$\mathbf{y}$, we define an error measurement using Eq.$\text{(3)}$$\eqref{eq:ref3}$ as follows

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{e}_{rr}& =\mathbf{y}-\hat{\mathbf{y}}\\ & =\mathbf{\text{y}}-A\hat{\mathbf{x}}.\end{array}\end{array}$$

Using the criterion from Eq. $\text{(4)}$$\eqref{eq:ref4}$ in matrix notation we have

$$\begin{array}{}\text{(6)}& \parallel e_{rr}{\parallel }^2=(\mathbf{\text{y}}-A\hat{\mathbf{x}}{)}^T(\mathbf{\text{y}}-A\hat{\mathbf{x}}),\end{array}$$

also, expanding the Eq. $\text{(6)}$$\eqref{eq:ref6}$ we obtain the following procedure

$$\begin{array}{}\text{(7)}& \begin{array}{rl}(\mathbf{\text{y}}-A\hat{\mathbf{x}}{)}^T(\mathbf{\text{y}}-A\hat{\mathbf{x}})& =\\ {\mathbf{\text{y}}}^T\mathbf{\text{y}}-(A\hat{\mathbf{x}}{)}^T\mathbf{\text{y}}-{\mathbf{\text{y}}}^T(A\hat{\mathbf{x}})+(A\hat{\mathbf{x}}{)}^T(A\hat{\mathbf{x}})& =\\ {\mathbf{\text{y}}}^T\mathbf{\text{y}}-{\hat{\mathbf{x}}}^T{A}^T\mathbf{\text{y}}-{\hat{\mathbf{x}}}^T{A}^T\mathbf{\text{y}}+{\hat{\mathbf{x}}}^T{A}^{T}A\hat{\mathbf{x}}& =\\ {\mathbf{\text{y}}}^T\mathbf{\text{y}}-2{\hat{\mathbf{x}}}^T{A}^T\mathbf{\text{y}}+{\hat{\mathbf{x}}}^T{A}^{T}A\hat{\mathbf{x}}.\end{array}\end{array}$$

To minimize the resulting expression in Eq. $\text{(7)}$$\eqref{eq:ref7}$ we perform

$$\begin{array}{}\text{(8)}& \underset{\hat{\mathbf{x}}}{min}\{{\mathbf{\text{y}}}^T\mathbf{\text{y}}-2{\hat{\mathbf{x}}}^T{A}^T\mathbf{\text{y}}+{\hat{\mathbf{x}}}^T{A}^{T}A\hat{\mathbf{x}}\},\end{array}$$

that is

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\frac{\partial }{\partial \hat{\mathbf{x}}}\{{\mathbf{\text{y}}}^T\mathbf{\text{y}}-2{\hat{\mathbf{x}}}^T{A}^T\mathbf{\text{y}}+{\hat{\mathbf{x}}}^T{A}^{T}A\hat{\mathbf{x}}\}& =0,\\ -2A^T\mathbf{y}+2A^{T}A\hat{\mathbf{x}}& =0,\\ A^{T}A\hat{\mathbf{x}}& =A^T\mathbf{y}\\ \hat{\mathbf{x}}& =(A^{T}A{)}^{-1}{A}^T\mathbf{y}\\ \hat{\mathbf{x}}& =A^{†}\mathbf{y}.\end{array}\end{array}$$

to obtain from Eq. $\text{(9)}$$\eqref{eq:ref9}$ the well know general expression of Linear Least Squares

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\hat{\mathbf{x}}& =(A^{T}A{)}^{-1}{A}^T\mathbf{y}\\ \hat{\mathbf{x}}& =A^{†}\mathbf{y},\end{array}\end{array}$$

where $\hat{\mathbf{x}}$$\mathbf{\hat x}$ is the estimated input variable and $A^{†}$$A^\dagger$ is the so-called Pseudo-Inverse matrix of $A$$A$.

Example 1: Worked Example with R

Consider the following model for data generation specified as

# Data generator functions
fdata <- function(t) {
  -2.50*t^2 + 1.75*t + 232.50 
}

we use $t$$t$ within range from $1$$1$ to $10$$10$ and generate simulated data using

xxxxxxxxxx
# Generate and plot observed data
x <- 1:10
ydata <- fdata(x)
plot(x, ydata, pch=19, panel.first = grid(),
     col='red', xlab = 't in seconds', 
     ylab = 'Observed data in Meters')
curve(fdata, from = 1, to = 10, col = 'darkred',
      lty=5, add = T)

The above code, results in the following plot which simulates some particle moving a certain amount of meter through the time.

given the observed output data $\mathbf{y}$$\mathbf{y}$, we wan to estimate the estimated input $\hat{\mathbf{x}}$$\mathbf{\hat x}$.

t	y
1	231.75
2	226
3	215.25
4	199.5
5	178.75
6	153
7	122.25
8	86.5
9	45.75
10	0

We suspect the model better fit is quadratic, that is

xxxxxxxxxx
# Fit Model
fmodel <- function(t) {
  c(t^2, t, 1)
}

Also we define support functions to compute pseudo inverse as

xxxxxxxxxx
# Pseudo Inverse compute function
pseudo.inv <- function(M){
  solve(t(M)%*%M)%*%t(M)
}

Now, we compute the measurement matrix $A$$A$ as

xxxxxxxxxx
# Coefficients estimation
x.hat <- pseudo.inv(A)%*%matrix(ydata)

to obtain the same values from the original model

xxxxxxxxxx
[1,]  -2.50
[2,]   1.75
[3,] 232.50

Under ideal conditions, if we observed a phenomenon, modeled it, and estimated its inputs, we would obtain its exact value, but this is not the case in reality. Consider now, the original model with additive noise as follows

xxxxxxxxxx
# Additive noise to original observed data
y.noise_obs <- ydata + rnorm(10,1,10)
A <- do.call(rbind, lapply(x, fmodel))

then to estimate the equation as in Eq. $\text{(3)}$$\eqref{eq:ref3}$ we perform

xxxxxxxxxx
# Coefficients estimation
x.hat <- pseudo.inv(A)%*%matrix(y.noise_obs)
y.hat <- function(xhat, t) xhat[1]*t^2 + xhat[2]*t + xhat[3]

where we use the following snippet to view the results

xxxxxxxxxx
# View results
curve(fdata, from = 1, to = 10, col = 'darkblue', lty=5, panel.first = grid(),xlab = 't in seconds', 
      ylab = 'Observed data in Meters')
curve(y.hat(x.hat, x), from = 1, to = 10, 
      col = 'darkred', add = T)
points(x, y.noise_obs, pch=19, col='blue')

Where the red line is the estimated model and the blue dashed line is the original and without noise model.

Conclusion

In many fields, such as reverse problems, identification systems for control problems, statistical learning models, and computer vision, least squares is a widely used method for estimating variables. This article provides a brief overview of the theoretical foundations and a working example. In our next article, we will examine the situation where it is necessary to estimate $\beta$$\beta$ coefficients in linear models for ordinary least squares regression problems.