Fitting a nonlinear-resistance model¶

Ohm's law says that voltage drop $d$ across a resistor is proportional to the current $i$ via $d = iR$, where $R$ is the resistance. Normally, we treat $R$ as a constant (the relationship is linear), but in fact this is not exactly true.

A real resistance actually changes as you increase the voltage drop $d$. (For example, this happens as the resistor heats up; eventually it melts!) One simple model of a nonlinear resistance is: $$ R = x_1 + x_2 d^2 . $$ Here, the resistance is constant plus a "small" correction that grows as we increase $d$. (We won't have any term proportional to $d$ if the current flows equally well in both directions, so that $R$ doesn't depend on the sign of $d$.)

One would normally get these coefficients $x_1$ and $x_2$ by experimental measurements. Such "fitting" processes lead directly into the next major topic in 18.06:

Exact fits¶

Suppose we put two voltage differences $d_1$ and $d_2$ across our resistor, and measure two resistances $R_1$ and $R_2$. This leads to the following $2 \times 2$ system of equations for the unknown coefficients $x_1$ and $x_2$:

$$ \begin{pmatrix} x_1 + x_2 d_1^2 \\ x_1 + x_2 d_2^2 \end{pmatrix} = \underbrace{\begin{pmatrix} 1 & d_1^2 \\ 1 & d_2^2 \end{pmatrix}}_A \underbrace{\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}}_x = \underbrace{\begin{pmatrix} R_1 \\ R_2 \end{pmatrix}}_b \; , $$

i.e. a system $Ax = b$. As long as $|d_1| \ne |d_2|$, the matrix $A$ is non-singular, and this has a unique solution: an interpolating polynomial going through the two data points:

In general, if you have $m$ points you can interpolate them in this way by a degree $m-1$ polynomial (called the Lagrange interpolating polynomial) in this way: you set up an $m \times m$ matrix system for the coefficients.

Note, however, that this is neither the most efficient nor the most accurate (given roundoff errors) method to compute an exact interpolating polynomial. There is a much better method called barycentric interpolation that is outside the scope of 18.06.

Inexact fits¶

Suppose we do more than just the bare minimum number of measurements, however. Suppose that we do $m > 2$ measurements. This gives the following system of equations:

$$ \underbrace{\begin{pmatrix} 1 & d_1^2 \\ 1 & d_2^2 \\ \vdots & \vdots \\ 1 & d_m^2 \end{pmatrix}}_A \underbrace{\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}}_x = \underbrace{\begin{pmatrix} R_1 \\ R_2 \\ \vdots \\ R_m \end{pmatrix}}_b \; , $$

Now, $A$ is $m \times 2$, and full column rank (assuming distinct voltages $|d_k|$), but of course it is not invertible for $m > 2$.

It will still have a solution (a unique solution!) if all of the $R_k$ measurements fall exactly on a quadratic curve $x_1 + x_2 d^2$, but in a real experiment there would be some measurement noise that would spoil this.

For example, let's suppose $x = (1,2)$ and we do $m=200$ measurements for $d \in [0,2]$, but that each measurement has a random uncertainty $\approx R \pm 0.1$.

The corresponding matrix $A$ is:

$A$ doesn't have an inverse, of course:

However, if we blindly do A \ b, it seems to give us quite a reasonable $x$, very close to the exact $x = (1,2)$ of the underlying data:

What is it doing? Let's plot the curve from the coefficients A \ b:

This $\hat{x}$ is not an exact solution: $A\hat{x} \ne b$: the "residual" (the "error") $b - A\hat{x}$ is not zero, as we can see by checking its norm $\Vert b - A\hat{x}\Vert$ above.

In the plot above, correspondingly, the black "fit" curve does not exactly match the data points. But it is pretty close!

Least-square fits¶

What A \ b is doing in Julia, for a non-square "tall" matrix $A$ as above, is computing a least-square fit that minimizes the sum of the square of the errors. This is an extremely important thing to do in many areas of linear algebra, statistics, engineering, science, finance, etcetera.

Above, we have a bunch of measurements $(d_k, R_k)$, and we are trying to fit it to a function $R(x,d) = x_1 + x_2 d^2$. There is no exact fit, so instead we minimize the sum of the squares of the errors:

$$ \min_{x \in \mathbb{R}^2} \sum_{k=1}^m (R_k - R(x, d_k))^2 \;. $$

Let's write this in terms of linear-algebra operations. We have our $m \times 2$ matrix $A$ above, and in terms of this we have $R(x, d_k) = (Ax)_k$: the "fit" function at $d_k$ is the $k$-th row of $Ax$. But then $R_k - R(x,d_k)$ is the $k$-th row of the residual vector

$$ \mbox{residual} = b - Ax $$

and we can see that our sum above is precisely the sum of the squares of the residual components, which is the square of the length of the residual:

$$ \sum_{k=1}^m (R_k - R(x, d_k))^2 = \Vert b - Ax \Vert^2 \; . $$

Here, $\Vert y \Vert = \sqrt{y^T y}$ is the length of a vector, also called the norm of the vector: the square root of the dot product with itself. (More specifically, this is the Euclidean norm, also called the $L^2$ norm. There are other ways to define a vector norm, but I think this is the only one we will use in 18.06.)

So, what we are really doing is minimizing the norm of the residual:

$$ \boxed{ \min_{x \in \mathbb{R}^2} \Vert b - Ax \Vert } $$

(Note that minimizing the norm and minimizing the squared norm will give the same solution $x$; which one we write is a matter of convenience.)

This is exactly what A \ b does in Julia (and Matlab) for a non-square matrix $A$, which is why it gives us a good fit above.

Least-squares fits and the normal equations¶

As derived in section 4.3 of the Strang textbook, minimizing $\Vert b - Ax\Vert$ or (equivalently) $\Vert b - Ax\Vert^2$ leads to the so-called "normal equations" for the minimizer $\hat{x}$:

$$ \boxed{ A^T A \hat{x} = A^T b } $$

These always have a solution. If $A$ is an $m \times n$ matrix, then $A^T A$ is $n \times n$, and:

• Usually in fitting problems, $A$ is a "tall" matrix with full column rank $n$, in which case $\operatorname{rank}(A^T A) = \operatorname{rank}(A) = n$ and $A^T A$ is invertible: the solution $\hat{x}$ exists and is unique.

• Even if $\operatorname{rank}(A) < n$, we still have a (non-unique) solution, because $A^T b \in C(A^T) = N(A)^\perp = N(A^T A)^\perp = C(A^T A)$.

Now (on the blackboard), I will show that we can derive the normal equations directly by 18.02. If we just take the partial derivatives of $f(x) = \Vert b - Ax \Vert^2$ and set them to zero, we find that $0 = \nabla_x f = 2A^T Ax - 2A^T b$ and the normal equations follow.

More polynomial fitting examples.¶

Suppose we are fitting m points $(a_k,b_k)$ to a degree-(n+1) polynomial of the form

$$ p(a) = x_1 + x_2 a + x_3 a^2 + \cdots + x_n a^{n-1} \; , $$

which leads to the $m \times n$ matrix

$$ A = \begin{pmatrix} 1 & a_1 & a_1^2 & \cdots a_1^{n-1} \\ 1 & a_2 & a_2^2 & \cdots a_2^{n-1} \\ 1 & a_3 & a_3^2 & \cdots a_3^{n-1} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & a_m & a_m^2 & \cdots a_m^{n-1} \\ \end{pmatrix} $$

and hence to the normal equations for the fit coefficients $\hat{x}$.

Let's generate 50 data points from a degree-3 polynomial $1 + 2a + 3a^2 + 4a^3$ plus noise, and see what happens as we change the fitting degree $n$.

Increasing the degree of the fit at first seems to improve things, up to about degree 6, but then it becomes more and more wiggly. It is still reducing the residual, but clearly it is "fitting noise" and the error actually becomes worse compared to the "real" underlying cubic model. This problem is called overfitting.

In between the fitted points, especially near the edges, the wiggles can actually diverge as we increase the degree, an effect related to what is known as a Runge phenomenon. Even if there is no noise, fitting equally spaced points to high-degree polynomials can lead to disaster if the underlying data is not exactly polynomial. Let's demonstrate this by fitting the smooth function $b(a) = 1/(1+25a^2)$ to polynomials at 50 points:

Be careful not to draw the wrong lesson from this. You have to be especially careful when fitting to high-degree polynomials, but this does not mean that you should never do it.

In particular, fitting smooth functions to high-degree polynomials can be a great thing to do, as long as you choose the correct points. (Separately, as mentioned above, the construction of the interpolating polynomial should technically not be done by this matrix method once you go to high degrees — you run into roundoff-error problems — but there are better methods like the barycentric formula.)

For example, here we fit the same $1/(1+25a^2)$ function as above to a degree-50 polynomial, with no Runge problems at all, by choosing the points to be Chebyshev nodes (which cluster together at the edges of the domain):

How this works goes far outside the bounds of 18.06, but is beautiful and fascinating mathematics. See, for example, this book and these video lectures by Nick Trefethen.

The great thing about fitting complicated functions to polynomials is that polynomials are usually easier to work with — finding roots, derivatives, and integrals of polynomials is easy, for example. A pioneering software package encapsulating this idea is chebfun, and a Julia implementation of similar ideas is ApproxFun.