Comment by wesm | Hacker News Reader

Comment by wesm | Hacker News Reader

Parent Full thread

If you spent a day doing some serious linear algebra in Python you might change your tune. * (multiplication) between NumPy arrays by default does element-wise multiplication (potentially with broadcasting), which is the desired default behavior. In R, for example, you can define custom infix operators so that

a * b

is elementwise but

a %* % b

is matrix multiplication. If you write down a complicated linear algebra expression, something like

A.T ! (B.T ! C ! B).T ! A

would be a lot friendlier to scientists than the current:

dot(A.T, dot(dot(B.T, dot(C, B)).T, A))

You might just say "well suck it up" but I've got to say that doing linear algebra in Matlab is a lot easier because the linear algebra that I do with pen and paper looks pretty much exactly the same as the corresponding code. On the other hand, Matlab is super clunky compared with NumPy at doing APL-style array processing with broadcasting operations, etc. In my work I tend to do more of the latter and less of the former but whenever I implement something with a lot of matrix multiplications it takes me a lot longer in Python to get things right.

Anyway, the point is: non-scientific Python folks need to take a walk in our shoes to gain an understanding of the challenges we face on a ongoing basis.

I'm having a hard time understanding your last statement. NumPy data types (dtypes) simply tell the ndarray how to interpret the block of data associated with it (the # of bytes per item, shape, and strides).

I actually agree with you. I much prefer operator overloading than unwieldy chained function calls. But it might be tough to convince the larger Python community to add it since it is only useful to a small part of the community.

The last comment about dtypes is that the size (shape in NumPy) should be part of the type; certainly element-type should be. A 1-D vector of double should not be the same Python type as 3-D array of characters or a 500x600 matrix. This creates havoc in a dynamic language. I once spent two days tracking down a bug caused by "*" multiplying two arrays instead of two matrices when I started using NumPy. Perhaps size is too much for a dynamic language to be part of type, but surely dimension and element-type should be reflected in Python type. It absolutely does not help that the documentation is littered with type-objects and object-objects; I appreciate this is how C implementations are, but for a beginner NumPy users, it is more than a little confusing.

I have a lot of respect for people who designed and built NumPy; for multi-dimensional arrays I don't see a better approach; but too much of influence of C implementation details seep into Python interface, types and operator overloading are only the beginning of the problems. I am having second thought about the suitability of dynamic languages for large scale, high-performance computing. Type annotation could help a lot but I see there is exactly zero interest in that for NumPy.

I see your point, and it is why I tend to like operator overloading in languages.

However, Python is a general purpose language, and scientific computing is a single domain. It is used in many different domains. Consider that there are many changes that individual communities would like, and if Python granted all of those requests, the language would be a mess. That's the challenge in designing a general purpose programming language.

> Python is a general purpose language, and scientific computing is a single domain.

I think this is shortsighted to the point of ignorance. "Scientific computing" here really means "performance-critical numerical computation on regular arrays". Basically all of the new things that people are doing with computers in the last five years and the next five years — machine learning and other statistics, software-defined radio, audio synthesis, real-time video processing, cool visual effects, speech recognition, machine vision, and 3-D rendering, and arguably Bitcoin — consist largely of performance-critical numerical computation on regular arrays. It's what GPUs are for. Five years ago, Numeric or NumPy was probably the best way to do that for a wide range of things, although a lot of people still use Matlab instead, and R deserves at least a mention. Today it's not clear. Five years from now there will be something much better than current NumPy, and it could be a better version of NumPy or it could be R or Matlab or Octave or something.

In short, "scientific computing" is not a single domain, but a set of capabilities increasingly important in many different domains.

I do research in high performance computing. I am well aware of the importance of performance critical, number crunching code that operates on large arrays. I wrote an entire dissertation based around better ways of expressing large computations on dense vectors and matrices for new multicores.

I am, however, self-aware enough to recognize that what I care about is a subset of what everyone in computing cares about. Scientific computing may be important in multiple places, but it is still a single domain. An important domain, sure. But Python is still used in many places were such concerns are not important. Arguments about the utility of number crunching on dense vectors and matrices are great. But the attitude I've seen in this thread is "the domain I care about is so important that my concerns should be elevated above the concerns of other domains." That is not going to fly when it comes to changing a general purpose language used in many domains.

I'm sure you know enormously more than I do about how to do dense-array number crunching. My claim, though, is that it's being applied much more widely than previously, not that there are new and better ways of doing it; so there are fewer and fewer places where such concerns are not important.

In the 1960s, there were those who claimed that the benefit of recursion was not worth the complexity it added to programming languages, except in certain special domains. In the 1970s and 1980s, we had the same argument about depth-first search — SNOBOL's and Prolog's backtracking feature. It turned out that the proponents of recursion were right, and the proponents of backtracking were probably wrong. (Although the backtracking feature of SNOBOL's child Icon directly inspired Python's generators, although implementationally they're maybe more similar to CLU iterators.)

Now, when it comes to matrix and vector manipulation, should we imagine it as a single domain, or as a feature that's useful across a wide range of domains?

I'm arguing for the latter.

And I'm arguing that's still one concern among many when it comes to a language used in so many different places. The arguments here are, I think, unsympathetic to the difficulties of designing and maintaing a general purpose programming language, in particular one with the design goals of Python.

If we define "general purpose" as "equally mediocre in all domains", then what you say makes sense.

If we define "general purpose" as "works well in as many domains as possible", then I'd argue that the correct response to an easily-remedied weakness in an important domain should be to address that weakness.

As you imply, the challenge here is to decide which requests to grant and which to deny, but the mere fact that a request is "domain-specific" (ignoring, for the moment, the questionable idea that _linear algebra_ is domain-specific) should not be enough to rule it out.

(As an aside, scientific / numeric Python is, I believe, one of its two or three most important application areas, and perhaps the most important historical factor in its success: that community has been championing Python since the days when people were using Perl for the web, or for anything else. Notice, for example, that Travis mentions working on SciPy in 1999---and the original "Numeric" package was written in 1995.)

Matrix multiplication isn't specific to scientific computing, it's a pretty common and large part of mathematics, and math is pretty general purpose.

I guess it would be possible to implement with a with-hack...

with numpy.doing_matrix:

    matrixy_multiply_goodness = m1 * m2

array_multiply_goodness = a1 * a2

operator overloading doesn't solve the problem because you frequently mix elementwise operations with matrix multiplication and division. Converting types between numpy arrays and matrices gets cumbersome

Reply on news.ycombinator.com