Least Squares in Python, JavaScript and MATLAB: Video Tutorial

This video tutorial first reviews the mathematical foundations of a special case of the Least Squares method: models that are linear in their unknown parameters. Then it implements the method from scratch in three programming languages: MATLAB, Python and JavaScript (using mathjs).

Mostapha Kalami Heris teaches the 39 minute 1 second video on the Yarpiz YouTube channel. The examples move from a straight line to models built from nonlinear features while keeping the unknown coefficients linear.

The problem

Assume that an unknown system takes an input vector $x \in \mathbb{R}^n$ and gives a scalar output $y$. Only a set of $K$ data points is known: the inputs $x_1, \dots, x_K$ and the matching outputs $y_1, \dots, y_K$. The goal is a function $y = f(x)$ that reproduces the data with the smallest possible error. The least squares method chooses the model that minimizes the sum of the squared errors between the model outputs and the measured outputs.

The key modeling choice is which features to include. Minimizing the training error finds coefficients for those features, but it does not establish that the model will predict new data well.

The linear model and its solution

The video collects the data into two matrices. Each column of $X$ is one input vector, so $X$ is an $n \times K$ matrix. The outputs form the $1 \times K$ row vector $Y$. The simplest model relates them through an unknown $1 \times n$ parameter vector $\theta$:

With measured data, $Y = \theta X$ is generally an approximation. Minimizing the squared residuals gives the normal equations. If the rows of $X$ are linearly independent, $X X^T$ is invertible and the solution can be written as:

The formula $X^{\dagger} = X^T (X X^T)^{-1}$ holds under that full row rank condition. The Moore-Penrose pseudoinverse is also defined for rank deficient matrices, where this inverse formula is unavailable. Using it gives a least squares solution with minimum coefficient norm when the minimizer is not unique. Prefer a numerical least squares method or pinv over explicitly forming the inverse, especially when features are nearly dependent.

Adding an intercept

A model of the form $y = \theta x$ always passes through the origin, so it cannot describe a line such as $y = x + 1$. The fix is to add a row of ones to $X$. The first parameter then acts as the intercept. With the data $x = (1, 2, 4)$ and $y = (2, 3, 5)$, the method finds $\theta = (1, 1)$, which is exactly $y = 1 + x$. When the last output is changed to 6, no line passes through all three points. The method then finds the best line, about $y = 0.5 + 1.36x$, and its outputs are close to the data.

Beyond straight lines

The same method fits any model that is a linear combination of known functions of the input:

Only the input matrix changes. Each column now holds the values $f_1(x_k), \dots, f_m(x_k)$ for one data point, which maps the inputs into an $m$ dimensional feature space. These functions are called regressors.

A polynomial of degree $p$ is a special case, with the regressors $1, x, x^2, \dots, x^p$:

The video fits a quadratic this way. Three data points with distinct input values determine one polynomial of degree at most two. With more points, least squares chooses the coefficients that minimize the squared residuals. The tutorial also replaces the quadratic feature with $e^x$ to show that the method is not limited to polynomials. Choose features from the problem and check their behavior on held-out data rather than assuming that more terms improve the model.

What the video covers

Theory (from 00:06)

The problem, the matrix form, the pseudo-inverse solution and the need for an intercept term.

MATLAB (from 06:53)

The video computes the pseudo-inverse with the formula and compares it with the built-in pinv function, which gives the same result. A script then fits a line to the data and plots the data points with the fitted line, using linspace for a smooth curve (from 09:35). After the general theory (from 13:36), a second script fits a quadratic and then a model with an exponential term (from 16:17).

Python (from 20:19)

The same two scripts are written with NumPy. The regressor matrix is built from np.ones and the inputs, the pseudo-inverse comes from np.linalg.pinv, and Matplotlib plots the data and the fitted curve. The results are the same as in MATLAB.

JavaScript (from 28:19)

JavaScript has no built-in matrix library, so the video installs mathjs with npm install mathjs and runs the code with Node.js (it can also run in a browser). The matrices are built with math.matrix and math.ones, and the pseudo-inverse is first computed step by step with math.transpose, math.inv and math.multiply. The video then shows a shorter way: math.evaluate accepts the whole formula in a MATLAB-like notation. The fitted model reproduces the outputs for the test inputs. A second script adds the quadratic term with math.dotPow.

Applications

Least squares is the basis of linear regression and curve fitting. It is used to fit models to measured data, to estimate the parameters of a system from its inputs and outputs, and to smooth noisy measurements.

Check the fit

Plot the residuals as well as the fitted curve. Systematic residual patterns can reveal a missing feature or an unsuitable model. Compare the scale of the features and check for duplicate or nearly dependent rows in the feature matrix. A low training error alone does not justify extrapolation, so test predictions on data that was not used to fit the coefficients.

Watch online

Play videoLoads from YouTube

Mostapha Kalami Heris, Least Squares in Python, JavaScript and MATLAB: Video Tutorial (YPTLS190514), Yarpiz, 2019. https://yarpiz.com/least-squares-in-python-javascript-and-matlab-yptls190514

@misc{yarpiz-yptls190514,
  author = {{Mostapha Kalami Heris}},
  title = {{Least Squares in Python, JavaScript and MATLAB: Video Tutorial}},
  howpublished = {Yarpiz},
  year = {2019},
  note = {Code YPTLS190514},
  url = {https://yarpiz.com/least-squares-in-python-javascript-and-matlab-yptls190514}
}