Can Your Static Type System Handle Linear Algebra? (2014) | Hacker News Reader

Can Your Static Type System Handle Linear Algebra? (2014)

yosefk.com

65 pointskushti65 comments

People keep talking about "you don't want to add meters to feet", but that's exactly the opposite of what you want to do. You want to add inches to yards or furlongs or flea wings, and you want to get it right.

What you don't want to do is add length to weight. The true type of your value is not inches, it's length, and expressing it in inches is a formatting issue.

Also, ideally, you don't want to add length to width, or Hertz to Becquerel (both have [1/sec] unit). This was my humble attempt to address this with a Scala Units library: https://github.com/adrianfr/scalau

If you want to find the periphery of something you DO want to add length to width...

No, you add the lengths of the sides. That's type-safe and generalizes from rectangles with edges parallel to the x and y axes to general polygons:

  perimeter = poly.edges.map(length).sum()

Height and width have (implicit) direction; length doesn't.

But yes, this will get hairy. For convenience, you want the ability to add heights, but that only makes sense in some cases. Two 100m high skyscrapers only make a 200m high skyscraper if you place one on top of the other.

Or maybe you want to find the height of an average skyscraper so you want to add and then divide. This is going too far maybe. At some point, one has to accept that not all logic bugs are type errors, and not all type errors are common enough in practice to worry about. I have no idea where that point is, but I'm sure it exists.

You might be interested in http://www.wolfram.com/mathematica/new-in-9/units/ - Wolfram's take on units.

When people say "you don't want to add meters to feet" it's a short-hand for "you don't want to add [a primitive data type representing] meters to [a primitive data type representing] feet [thereby getting a nonsensical result]." Whether "inches" is the "true type" of my value depends on the context - what invariants am I trying to enforce by giving it that type?

Don't skip the comments. I think there is more there than in the article, honestly.

Of course there is. What there isn’t, as far as I can see, is a working solution.

There is, more or less. Mine, in fact. But it's ugly. It's repetitive. If someone gave me that solution, I would think "this is way too unreadable to ever verify its correctness". My mistake was using C++98. I think C++11 would handle this better, but I never got around to trying.

There's no reason type(X'X) should be fixed. If you allow yourself to consider the case where each entry in the matrix has different type everything is fine. You just end up with matrix product rules like:
[type(a) type(b)][type(c); type(d)] = type(a)type(c) + type(b)type(d)

Which puts constraints on the types of the various entries in the product, but doesn't restrict them to being the same in each matrix or vector. This is totally reasonable - in Ax + b = y I'm restricted to type(A)type(x) = type(b) = type(y) for scalars. Why wouldn't we have a similar generalization of constraints in the matrix case? The model that arrays of objects must all be the the same type simply doesn't fit with the analysis in this case. The author suggests this:

> But I just can't imagine the monstrous types of products of the various elements ever eventually cancelling out as nicely as they apparently should – not in a real flesh-and-blood programming language.

Well that's how it should work. You just have to imagine it :)

Perhaps the sanest approach is we type every column of X and Y separately...

Typing every column of a matrix is not sufficiently general. Each element of the matrix can have its own type.

Suppose you're solving a linear algebra problem in two unknowns. Say it's an electromechanical system, and we want to find a voltage and a pressure that gives us some desired flow rates. So let x [=] (V, Pa), for example, and let y [=] (mA, kg/s).

Then the matrix-vector product Ax = y requires the units:

    [ mA/V   , mA/Pa   ] * ( V  ) = ( mA   )
    [ kg/V/s , kg/Pa/s ]   ( Pa )   ( kg/s )

Of course those units in A have nicer reductions, but it doesn't teach us anything to figure them out. The point is that the matrix A represents the system design, so the units don't have to have familiar names, let alone be consistent within a column. If the system involves Johnny reading a value off a meter, and then phoning somebody up to tell them to pedal a stationary bike a little faster, that process has some kind of units associated with it, and the matrix A must include those units in order to represent the system.

Typing every column of a matrix is not sufficiently general. Each element of the matrix can have its own type.
While this is true, letting every matrix element have an arbitrary type is too general. The matrix' determinant should have a well defined type to be invertable. In your example mA/Pa * kg/V/s and mA/V * kg/V/s are compatible, so it works out.

It's not true that each element can have its own units; given an MxN matrix there are only M+N-1 degrees of freedom (basically you can choose arbitrary units for each element in the first row and then a single multiplicative factor between the first row and each other row). Otherwise you'll find that either left or right multiplication won't have any compatible vector types.

I think this is the strongest argument for dynamic typing I have ever read.

Wow, if that is so, then I am completely sold on static typing. I do not get how dynamic typing gets into the picture.
Having a type system that can catch some errors at compile time but not all is better than having a type system that doesn't catch any.
Having a static type system doesn't mean that you can't use dynamic typing as much as you want (e.g. in Scala you can make all variables of type `Any` and use type refinements whenever you want to call a method).

You seem to be under the impression that I'm a proponent of using dynamic typing to the exclusion of static typing. That's not the case (I'm an ObjC fan). I just hadn't seen a situation before where the argument for one vs the other in a particular situation hinged on something so fundamental rather than on a design artifact.
EDIT: Removed comment about scala because I was confusing it with clojure. My bad.

That's the thing. There is no argument for one vs the other.
There is a statement of the following fact: "It is hard to express heterogeneous array type in most statically typed programming languages and verify operations at compile time." Well, duh. AFAIK this would require dependent types* (not something that you get for granted in mainstream languages).
But can we check types in runtime? Yes. Does it mean that it can only be done in a dynamic language? No. It can be easily done in any language that has integers and arrays. Even x86 assembly is enough.
* Scala's type system should be enough, if we use shapeless.HList's for column and row types and then a ton of implicits to resolve the resulting type and check that everything is sound.

What's striking about this is that this constraint is essentially a type constraint. To the degree that existing type systems (and the checkers therefore) make this ugly, difficult, or impossible to express that surely points to an opportunity for improvement.
Note that, contrary to some on this page, I read this as speaking entirely about computation whose structure is known at compile time. Obviously, if the structure is based on run-time information then you probably need dependent types to guarantee statically that nothing can go wrong.

Are there dynamically-typed systems that actually do this? That check the units of matrix and vector elements at runtime without checking every single time you want to multiply two numbers together?

I imagine most languages with a matrix library uses opaque types implemented in C or similarly fast language. So just one typecheck would need to be done.

I'm sure they check that the matrix doesn't contain a string, or that you're not mixing 32-bit and 64-bit floats, but checking user-defined units as in the article, where each row or column may have a different one seems almost as hard in a dynamically-typed language, unless you do the check at the very bottom level. And that would be unimaginably slow.

Oh ok, yeah I refered to 32/64bit floats / ints. I didn't like the notation in the article, so I didn't bother understanding what he was talking about.

I'm a statistician by training, though these days I do a variety of work as an engineer.
I'm a huge proponent of static typing, but I have to say that there's a reason that the predominant languages used by statisticians and scientists are generally dynamically-typed (Python, R, MATLAB, etc.) - and it's not that these people simply "don't know how to deal with static types"[0].
I haven't used Haskell in a few years, so it's a bit rusty, but I'll use it as an example, since it's the most widely-known strict static type system. Hindley-Milner type systems like Haskell often don't allow dependent types[1], which means that a 2D vector and 3D vector are the same type unless their dimensions are known at compile-time (which is oftentimes not the case). This is a big problem, because I may want to run a regression on an unknown number of variables, but also constrain the number of variables in alternative regressions to use the same number of variables.
In Haskell, you might solve this problem with a tuple, except then (1) your vector elements may be of the wrong type (bad!), and (2) You can only have vectors of length 63 or less.[2]
By contrast, in MATLAB, it's trivial to declare a number that may also be a vector, but act like a scalar for certain operations, but also be correctly "incompatible" with vectors of incorrect dimensions[3].
Yeah, you can do this in Haskell/Scala/Go/<insert-language-here>. But it's really, really annoying.
There is one assignment that I had in school which was done in MATLAB, and just for fun, I tried reimplementing it in a few other languages (Haskell, Python, Go, and a few others). I ended up giving up in frustration with all three of these, whereas the original MATLAB assignment only took a few hours (most of which was figuring out the math behind it, not writing the code). My Python attempt was the one which came the closest, but even then it was just so hideous that I couldn't bear to work on it any longer[4].
[0] I haven't used Julia, so I can't say whether the optional typing there solves this problem or not.
[1] At least, Haskell didn't as of a couple of years ago, which was the last time I checked. Though it seems this is still an open problem: https://wiki.haskell.org/Dependent_type#Dependent_types_in_H...
[2] https://downloads.haskell.org/~ghc/7.0.4/docs/html/libraries...
[3] If you think that this sounds easy in your language, try generalizing this where `v` can be an element in any vector space, not just the sorts of vector spaces that we're used to dealing with as computer scientists.
[4] Yes, I am aware of the irony in saying saying that this Python code was more hideous than the equivalent MATLAB code. MATLAB is awful in so many ways, but at the end of the day, it's a domain-specific language, and it tackles its use case in a way that no general-purpose programming language I've seen has.

In response to [0], yes Haskell still lacks dependent types. Fortunately there is research going on. Liquid Haskell allows you to use refinement types to help further constrain things. Anyone who is interested in full dependent typing should look at Idris.
http://www.idris-lang.org/

> I'll use Haskell as an example, since it's the most widely-known strict static type system. Hindley-Milner type systems like Haskell often don't allow dependent types, which means that a 2D vector and 3D vector are the same type unless their dimensions are known at compile-time (which is oftentimes not the case).
Implementation of dependent types in Haskell, with the specific example of dimension-parameterized vectors, is addressed in [0].
[0] https://www.fpcomplete.com/user/konn/prove-your-haskell-for-...

What was the problem in [4]? I'll agree that there are a couple nice things about the core matlab syntax, but things tend to get very terrible very quickly once you start to actually do anything with it.

Now I'm curious: what was this problem? I can't imagine any program in Matlab that one can't simply transcribe into Python almost mechanically.

Here is some stuff written on Julia: https://github.com/JeffBezanson/phdthesis/commit/9832392cce6...

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) 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.

F# units of measure to the rescue:
https://msdn.microsoft.com/en-us/library/dd233243.aspx
let convertCtoF ( temp : float<degC> ) = 9.0<degF> / 5.0<degC> * temp + 32.0<degF>

F# units can certainly handle the problem given a fixed matrix size, but for a general method handling arbitrary dimensions you'd need variadic measure arguments, which aren't supported (and would probably require quite a bit of bookkeeping for this particular problem, in any language).

Look at the "Eigen" [1] library for C++. It's totally possible to do linear algebra in a staticly typed environment. You just need to use generic types. However it isn't as pleasant as doing it in a dynamicly typed environment like MATLAB.
[1] http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra....

as far as I know, Eigen is type-safe only with regards to dimensions of vector spaces. It can check at compile time that a dot product is taken over two compatibly-sized vectors, or that a matrix-matrix multiply makes sense dimensionally.
It doesn't have a notion of vector spaces with units. It can't keep track (not Eigen, at least) that a vector [m, kg] cannot be inner-product'd with itself unless its through a norm that produces the correct dimension, otherwise the units are m^2 + kg^2.

Isn't that the default behaviour for most linear algebra systems? (e.g. MATLAB). Maybe I'm missing something but I don't really care about compile time checking and it seems an impossible task anyway with variable/dynamic sized inputs.

Yes, it is impossible for a compiler to know the size of a matrix whose size is determined at runtime. But a lot of times you do know the size of matrices and not only is it more type-safe to have the type of a matrix be parametrized by its size, it also allows for more efficient code generation. If you don't know the size of a matrix and vector, you necessarily need a loop. If you know your matrix is 3 by 2, you may just, at compile time, unroll the loop completely.
I don't know which default behavior you mean.

It's impossible for a compiler to know the size of a matrix whose size is determined at runtime, but it's not impossible for a compiler to prove that whatever the size your logic is still correct (or at least doesn't break any of the invariants you've established). Depending on the nature of those invariants and the complexity of your logic, you may need dependent types.

I can't answer the question directly, because as a web engineer I haven't once in my career actually had to do linear algebra in my code. That said, the entire mindset behind this article strikes me as fundamentally wrong-headed. The mindset is “either you have static types, or you don't.” From that perspective, if you can think of any situation in which you can't encode the problem into static types, then you might as well throw the whole idea out the window and just hope things work at run-time. In the real world, this isn't at all how it works.
The last couple of decades have been a story of type systems slowly but surely expanding the class of problems for which type-safety is practical and worthwhile. There are still problems for which this isn't the case, and there will likely still be problems like this for a long time. Proponents of static typing aren't arguing against this fact; they're arguing that for all the myriad cases where static types can help, we should use them. For the cases we can't figure out, fine, just use ints. But for those functions, write a lot of tests, and be sure to rewrap those ints in semantic types as soon as possible.
In short, whether or not my static type system can handle linear algebra is completely irrelevant to me. To be worth using, my type system doesn't have to solve every problem. It only has to solve enough problems to provide more benefit than cost. In my experience, it does that by multiple orders of magnitude.

> solve enough problems to provide more benefit than cost
What kind of problems do you have in mind where a statically typed lang would provide more benefit than cost over modern dynamically typed langs (erlang, clojure, etc.) or just in general?

In my experience, most application-level domain problems benefit significantly from strong static types. My most salient experience in the area is from working on a data analysis pipeline web app written in Haskell. The type checker regularly caught bugs that would have taken me hours to deal with if I had had to track them down after something broke at runtime.
On the other side of the equation, again in my experience, the actual cost of a good static type system is minimal. Type inference means that you don't waste time on annotations, and the time you spend writing up the types themselves is time you would be spending anyway, since the primary difficulty is understanding the domain and the shape of your data.

Thanks. It seems you are more concerned more about the software development experience and in that regard you see benefits. But I'm more curious and concerned about the benefits in solving real world problems. For example, in what way would Haskell be more beneficial in solving the problems that Erlang currently solves at Ericsson.

I generally feel that solving problems with software is the easy part.
Solving in the way that is extensible, maintainable, bug-free as possible, and efficiently is the important part.

> It seems you are more concerned more about the software development experience and in that regard you see benefits. But I'm more curious and concerned about the benefits in solving real world problems.
I don't see the distinction.

The software development experience is the real world problem. Programming languages are force multipliers. My argument is that statically typed programming languages are force multipliers with bigger coefficients.
I doubt that Haskell would be more beneficial for Ericsson, because Erlang was literally designed for their exact use case. I would be more inclined to pitch Haskell to the >99.9% of companies who don't have battle-tested custom programming languages.

Replace Ericcson with WhatsApp. What benefits would WhatsApps product see if they had used Haskell?
The software development experience is exactly that, an experience. You see problems in it, but others do not in regards to static and dynamic typing.

Haskell has more options WRT concurrency and parallelism than Erlang does at a single-machine level and the options for clusters are widening quickly.
Examples:
https://github.com/bitemyapp/blacktip (simple MVar based synchronization producing a much faster program than the Erlang version with zero optimization effort)
http://haskell-distributed.github.io/
http://hackage.haskell.org/package/courier
But don't take my word for it. I'd rather train coworkers eager to learn and retain an institutional competitive advantage than convince HN users :)

The software development experience is not only about subjective opinions ("I like static typing better!") but also about reliability and bug-catching.
Forget Haskell. I'm sure you know this is flamewar material, but the purported advantages of statically typed languages over dynamically typed languages are that it's easier (less bug prone) to refactor code written in statically typed languages, and also that there are fewer tests you have to manually write, since the compiler "writes" those tests for you, so to speak, helping write correct software with less effort. In general, the argument goes "all else being equal [1], it's easier to write correct software using a statically typed language". As an end-user, the software development experience affects you even if you are not a programmer!
I'm sure you won't agree with this. You could as well ask Paul Graham "what benefits would WhatsApp see if they had used Lisp?", and his answer would likely be as unconvincing to you as mine in this post.
----
[1] availability of libs, performance, and not competing with a custom language designed specifically for the problem at hand.

Your question was already answered: "regularly catching bugs that would have taken hours to deal with". That is the benefit nbouscal believes WhatsApp or pretty much anybody else would see. You can argue whether sophisticated static type systems accomplish such a thing, but it seems inarguable that it would be widely beneficial if they do.
I'll add my own (related) real-world benefit that I feel when I'm using a good static type system: I find it easier to write good unit tests when the types of data being passed across units has been statically limited.

The problem is that linear algebra is probably one of the most basic things you can do with a computer. If static typing isn't capable of offering useful help or even just has a few weird wrinkles, then something is very odd.
For what it's worth, I've seen LA done in C++ with various templates and libraries, and so it's certainly possible in various ways.

It wouldn't surprise me at all if static typing was capable of offering useful help. If it can't today, I wouldn't be surprised to see us make progress over the next few years. I think interesting conversations could be had about why this problem is harder to model using type systems than other problems are. Where I get confused is when this is used as an argument against using static typing in general. By analogy:
“Hammering nails is probably one of the most basic things you can do as a carpenter. If lathes aren't capable of offering useful help, then something is very odd.”
I don't mean to mock, but it doesn't seem odd to me at all that not all tools are good for all purposes. What seems odd to me is trying to argue against the use of a tool that has proven itself in many different situations by pointing out that there exist situations it isn't useful in. Sure, that's true, but so what?

but a programming language isn't a hammer or a lathe. its a workshop that is used to produce hammers and lathes and every other tool you will need.
the criteria for judging the workshop is rather different than the criteria for judging the lathe.

PL's as tools is a half-accurate analogy. If you're going to stretch it that far, shouldn't you use something more accurate? Or better, not use an analogy for something that detailed?
Say, PL's are materials. They are something that useful things "tools" are fashioned out of. The properties of a material affect the entire thing built out of it, in modest but pervasive ways. Some programs are built in more than one language, this could be seen as an object built with different materials for different parts. Now, whether mainstream languages are analogous to different iron alloys, wood, concrete or plastic is a different discussion. (And maybe another exercise in taking an analogy too far...)
Not sure where I saw this analogy first, since I can't find it on LtU...

> I think interesting conversations could be had about why this problem is harder to model using type systems than other problems are.
We can summarize this in ten seconds: More equations hold between combinations of units than combinations of types constructed using type constructors. For example, m kg m^-1 = kg, but no such relationship exists between (A, (B, A -> R)) and B. Proving these sorts of equivalences adds complexity to type systems.

I think you glossed over the point of the article in general. The author isn't asking whether Linear Algebra can be encoded into a static type system, hell, we had Fortran over 50 years ago! What the author is asking is whether the static type system can reasonably check the units of each member within a matrix or vector. Let's go over a short example:
Assume that you have a matrix A, which contains several rows corresponding to measurements, and several columns corresponding to the partials in x (such as in the generic parametric least-squares problem => A * x - b = r, where we want to minimize r^2). What the author is asking is whether or not the static type-checker can manage units and conversions within each element of the matrix. Let's say the first 10 rows of A correspond to units of metres, while the next 10 correspond to angular units (i.e. radians or degrees). The solution of such a system is:
x = -(A' * A)^-1 * A' * b
What the author desires is a way to enforce that the units of the elements of A are maintained, similar to if you had:
LengthInMetres a = LengthInMetres 1 + LengthInInches2;
Where addition between LengthInMetres and LengthInInches would perform the appropriate conversion if need be, or else you would raise an error for a type contract violation. Currently, there isn't any system that encodes the types of each row, column, or even the individual elements within a matrix, which is why the author is raising the question whether or not it's possible for a type system such as Haskell's to do so.
Right now, the best I've seen is linear algebra packages that can template over the internal type, but represents everything as doubles, or a number type of some kind. Numpy has some ability to modify internal data types (such as representing a matrix where each row is a tuple type or some other), but there is:
1. No static enforcement of type contracts in Numpy (more specifically Python)
2. Numpy does not provide the full set of matrix functionality if types are not the standard numeric types (obviously, you can't take the determinant of some arbitrary `object`, even if internally it represents an integer in metres).
Overall, what the author raised is an interesting question, and makes me wonder if we'll ever get to such a point in the future. Of course, seeing as most matrix libraries are built on top of BLAS / ATLAS / LAPACK / MKL, it's doubtful that we'll see any high-performance linear algebra with smart type systems.
EDIT: Formatting

The type in actual numerical implementations tends to be 'double' or 'real'.

>>I can't answer the question directly, because as a web engineer I haven't once in my career actually had to do linear algebra in my code.
Wow, really? I am a web developer and I have to do tons of math. Not just the WebGL/VR stuff, but the CRUD stuff too. Reports and queries and, yes, unit translations.
Like, what do you do? What are you making that doesn't need math? I've never not needed math. I wish I knew MORE math.

He did not say "math", he said "linear algebra", which is much more specific and specialized.
When was the last time your "CRUD stuff" involved solving a set of linear equations?

No, no, I think it's one of the most generally useful fields of applied math. You need it for best-fit functions, at the very least.

It's incredibly rare to not use a library for this kind of math unless you actually work on mathematical software.

Look, the person I was replying to clearly meant to express "web engineer, therefore not mathematician". That's my only point.
And it does you no good to have a linear algebra library if you don't know how linear algebra works.

Given that the person you were replying to said, "I apply category theory and type theory all the time, and I'm sure the time I've put into studying algebraic topology and abstract algebra has shaped how I think about problems[,]" I don't think it's clear at all that they meant to express "web engineer, therefore not mathematician."

That doesn't really mean anything

We can pick at the semantics, but they don't really matter. My point was that the pragmatic function of the above is quite clearly at odds with an attempt to distance the poster from mathematics.

I didn't say I don't use and need math: I apply category theory and type theory all the time, and I'm sure the time I've put into studying algebraic topology and abstract algebra has shaped how I think about problems. I actually really enjoy linear algebra as well, and always enjoy the opportunity to learn more of it.
What I said was that I've never had to do linear algebra in my code. Most of my work has been building web apps. Most of the domains that I've built those web apps for don't involve mathematical computations more complicated than multiplication, and when they have, I've usually just passed the task off to pre-existing libraries.

Reply on news.ycombinator.com