This isn't really a problem with the type system. Basically, the author figured out that the normal equations with Vandermonde matrices (https://en.wikipedia.o

This isn't really a problem with the type system. Basically, the author figured out that the normal equations with Vandermonde matrices (https://en.wikipedia.org/wiki/Vandermonde_matrix https://en.wikipedia.org/wiki/Vandermonde_matrix) are horribly ill-conditioned and a terrible tool for interpolation. Though, I actually like the authors example because it's a nice way to give physical units to what sometimes causes ill-conditioning. But, in short, I don't think there's anything wrong or complicated about a type system tracking units in this case. Vandermonde matrices will absolutely have radically different units in their elements, which is a huge disincentive to use them. Now, what's a better way to interpolate a series of points? Basically, use a better polynomial basis. In the authors example, we use the polynomial basis [x 1]. Or, in the higher order case, [x^4 x^3 x^2 x 1], etc. A better way to do this is to use a polynomial basis that's more stable, such as Legendre, Chebyshev, or Bernstein polynomials. They all have different properties, but if we look at something like the second order Bernstein polynomials, we have [(1-x)^2 2x(1-x) x^2], we notice that each polynomial has the same units. Then, if we do form the normal equations, we don't have such a huge disparity in the units of the elements. In addition, forming the normal equations are a terrible way to solve these problems. Basically, if we want to find a least squares solution, we're much better off using something like a QR factorization because things like Givens rotations are much more numerically stable. Alternatively, we can use iterative methods like LSQR or LSMR, which also compute the least squares solution without explicitly forming the normal equations and thus avoiding the condition number problems. Finally, even if someone wants to be neurotic about units, which is generally good practice, but a pain to do every time, we're almost always better off doing dimension analysis at the operator level and not the element level. For example, if I was solving the heat equation, I might form a matrix with the finite different operator or a stiffness matrix. However, I would know that the stiffness matrix, in this context, would have units of L(m^2/s,K/s) where L means that we have a linear operator from something with units of m^2/s to K/s. Technically, we can work backwards from these two to see that our elements have elements of K/m^2, but that's not as useful as the operator formulation. And, I claim this because if we write things out in the operator form, it makes it easier for us to figure out the units of the adjoint. In the simple case with the L^2 norm, we know that the adjoint has units of L(K/s,m^2/s). However, if the inner product we use has units, this answer will absolutely be different. Personally, I find that easier to track when I write things out like this. Anyway, I do think dimension/unit analysis is important, but static type systems are sufficient to track it and it's really not too bad to do so.