Thinking in an array language
github.com
github.com
Story time: 2 years ago, right before 1st lockdown, I landed a client. A data scientist who had already implemented an algorithm which dealt with matrices, implemented by a previous programmer he hired. Said implementation was in Python, no more than 10 line altogether, which performed well when matrix size where small, like 10x10. But the problem was, his real work need it to have matrices of size 10^6 x 10^6. Not only the Python implementation had to be ran on a beast of server with 4TB of memory, it also took 3 days to finish. And while the algorithm was small in size in Python, its explaining paper was 4 pages in total, which took me 1 week to understand. And then an entire month to implement in C. But in the end, when all was said and done, the run time when used with real data was only 20 minutes and consumed only 8 GB of memory, though it did required at least 16 virtual processors.
Hence my question, in the end performance is what it matters, not number of lines.
a:100 0N#?10000 / 100x100 double matrix
\t:1000 a(+/*)\:a / Total of 1k runs in ms
1514
#include <stdio.h>
#include <time.h>
size_t monoclock(void) {
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC, &ts);
return 1000000000*ts.tv_sec + ts.tv_nsec;
}
int main() {
const size_t n = 100;
double a[n][n], r[n][n];
size_t t = monoclock();
for (size_t iter=0; iter<1000; iter++) {
for (size_t i=0; i<n; i++) {
for (size_t j=0; j<n; j++) {
double sum = 0;
for (size_t k=0; k<n; k++) sum += a[i][k]*a[k][j];
r[i][j] = sum;
}
}
}
printf("%lld\n", (monoclock() - t)/1000000);
}
Also the scientist was leading a team, so my program would've been used by at least 20 people, 10 times per day. That was why Python one was a no go for them from beginning.
The C version wasn't GPU-targeted though from your description. I'm curious what other implementations would be capable of, for instance julia, maybe gpu-targeted.
It would be better to do something with the computed result so the compiler does not remove the computation, e.g. print a randomly selected value.
I also benchmarked the fixed C code against Python and "Python" (which uses OpenBLAS under the hood in my case) was 10 times faster:
import time
import numpy as np
a = np.random.rand(100, 100)
t = time.perf_counter()
for _ in range(1000): a @ a
print(time.perf_counter() - t, "seconds")
a =: ?100 100$0
<.1e6 * 1000 (6!:2) '+/ .*~ a'
269
a =: ?200 200$0
<.1e6 * 1000 (6!:2) '+/ .*~ a'
1796Even the older Mathematica language blows both K and C out of the water for readability and performance:
m = Table[RandomReal[], 100, 100];
t = RepeatedTiming[MatrixPower[m, 2]];
First[t]*1000*1000
18.0245
Oh, and Mathematica's matrix power function is:
- Parallel! If I use a 10K x 10K matrix as an input, it uses about 75% CPU on my 8-core laptop. It can complete a single multiplication in 5.3 seconds. For laughs, try that with either K or C and see what you get...
- Extends to negative or fraction powers.
- Has optimisations for applying the matrix power directly to a vector.
- Is extensively documented, unlike the terse K snippet or the hand-rolled C code: https://reference.wolfram.com/language/ref/MatrixPower.html....
Essentially what I'm trying to say is that terseness is an anti-pattern, and doesn't even begin to approach the utility of a well designed high-level language intended for teams of collaborating humans.
% python matrix.py
Timing 10 squares of a random 10000 x 10000 matrix
97.3976636590669 seconds
python matrix.py 364.41s user 8.10s system 379% cpu 1:38.25 total
% julia matrix.jl
Timing 10 squares of a random 10000 x 10000 matrix
97.787679 seconds (31 allocations: 7.451 GiB, 0.33% gc time)
julia matrix.jl 405.34s user 8.13s system 375% cpu 1:50.09 total
These folks are really onto something, but I think they get sidetracked in the (admittedly very very fun) minutia of the languages and lose sight of the crucial insight in re: mathematical notation, to wit: it's a means of human communication.
For APL or K to get out of their niches would require, I am convinced, something like a tome of documentation of a ratio of about 1.5 paragraphs per line of code. That would give us mere mortals a fighting chance at grokking these tools.
A similar problem plagues the higher-order stuff they're pursuing over in Haskell land. I know e.g. "Functional programming with bananas, lenses, envelopes and barbed wire" and "Compiling to Categories" are really important and useful, but I can't actually use them unless some brave Prometheus scales Olympus and returns with the fire.
Stuff dribbles out eventually. Type inference and checking have finally made it into the mainstream after how many decades?
You can download a free PDF copy of Mastering Dyalog APL by Bernard Legrand which is 700+ pages, from here:
There was a better example a couple of weeks ago here in a thread, someone had done a bit of APL or K for Advent of Code or something and posted a line, and someone else broke it down and explained how it worked. I spent an hour just now with Algolia trying to find it but I failed. :(
Holm's Klong docs have a good example in https://t3x.org/klong/klong-intro.txt.html where he explains how a table formatter function works.
Mathematical equations are usually embedded in papers that explain them. (I mean, I've read papers that were basically equations one-after-another with just scraps of interstitial prose, but they were heavy going.)
It wasn't a couple of weeks ago, but I did that for a line here: https://news.ycombinator.com/item?id=30463080
Or could it have been on an Advent of Code link? There have been some explanations in the answers mega-threads on Reddit. Anyway, yes I agree more explanations would be benficial - and I think there would be room for an animated explainer website with small blocks representing the array elements, coloured by how they are grouped by each primitive operation, and visually showing them moving round and splitting and combining. Such a thing would make a lot more sense for an array language than for many languages.
Your explanation was fantastic, and yeah I think the availability of more information like that would go a long way towards lowering the barrier for people to pick up these languages. Even if they don't actually use APL and it's ilk they can still get a better idea of how to use things like Numpy by getting familiar with the concepts and idioms that support them.
(I forgot to mention "A History of APL in 50 Functions" https://www.jsoftware.com/papers/50/ I just started working my way through that and it's been very helpful.)
> and I think there would be room for an animated explainer website with small blocks representing the array elements, coloured by how they are grouped by each primitive operation, and visually showing them moving round and splitting and combining. Such a thing would make a lot more sense for an array language than for many languages.
That would be pretty awesome. Reminds me of Guo's Python Tutor https://pythontutor.com/
https://docs.google.com/document/d/1W83ME5JecI2hd5hAUqQ1BVF3...
https://github.com/tlack/b-decoded
https://chat.stackexchange.com/rooms/90748/conversation/ngn-...
They're not 1.5 paragraphs per line, but enough to give a taste of the implementation style.
matmul: (+/*)\:https://aplwiki.com/wiki/Iverson_notation
https://aplwiki.com/wiki/Comparison_with_traditional_mathema...
If they wrote in an array language instead of LaTeX, not only would it make writing papers easier (+/ is shorter than either \Sigma or \sum), but it would be trivially reproducible, due to being an executable notation.
Why is that?
For example 'a*b' is any number of 'a's followed by 'b'.
How else would you concisely state that?
Presumably people who hate array languages think all 3 character regexes should instead be big nested loops, so they are "readable".
When writing a program, I prefer to use a PEG, giving the less compact notation `'a'* 'b'` but also letting me say `'a'* b` and define b as its own rule, including recursion for the useful cases. It helps that it's more powerful, being little more than a formalization of the post-regular strategies used in Perl-style 'regular' expressions while embracing recursion.
For '/' in vim, grep, wherever? Yeah regex is fine, that's what it was designed for.
"a"* "b"
Could you explain this? I don't quite understand what the problem is. Do you mean something like sed's regex substitute command?
If you use syntax like "a"* "b" then it's really obvious - the stuff in quotes is literal text, everything else is control.
Lots of formats make the same mistake, e.g. YAML.
Anyway I'm not sure the same is true for K. At least for the given example the for loop was not exactly difficult to write.
Everything is unreadable until you learn how to read it.
The moment you need to run complex group bys and store non-contiguous data efficiently, it gets awkward pretty quick.
On the other hand, operations on dense data is pretty cumbersome in SQL. You can only do so much with limited support of proper scan algorithms, merge joins or bolted on window functions.
Please somebody combine APL with SQL and you win the programming language wars.
There's a lot of database software built on array languages, with kdb and Jd being the most prominent as far as I know.
[0] https://aplwiki.com/wiki/Key
[1] https://code.jsoftware.com/wiki/Vocabulary/slashdot#dyadic
kdb+ and q fit this description. Docs here: https://code.kx.com/q/basics/qsql/
Here's an example.
In SQL:
SELECT stock, SUM(amount) AS total FROM trade GROUP BY stock
q)select total:sum amt by stock from tradeI was fortunate to have a teacher whose database course drilled relational algebra into us. This was in the 70s, shortly after Codd's paper, and well before SQL was invented, much less established. Now I think about much computation algebraically (and often functionally). But I do see that this is "unnatural" for many, including students, having taught databases for several years.
SQL reflects this. I often see students writing nested subqueries, because that is more procedural, where joins would be a cleaner choice. A colleague of mine wrote a paper many years ago, pointing out that thinking procedurally is more "natural" for many: https://dl.acm.org/doi/10.1145/319628.319656. But thinking set-at-a-time instead of one-at-a-time is a valuable skill, not that far off from thinking functionally.
- Data modeling - Relational algebra - SQL - DBMS architecture (buffering, btrees, query processing, index selection, ...) - Query optimization - Transactions
A few assignments have students write an in-memory relational algebra implementation, and then use it to write queries. A typical query is a deeply nested set of function calls (as a "one-liner"). And only then to we get to SQL. The hope is that RA is so ingrained, that the connections to SQL are easier to see. And this background is also really useful in understanding query optimization.
All of this material, including assignments, is available online (this is my own webserver): http://geophile.com/115.
Give me 'any join' operation - ie. just select the first value instead of all, and I'll happily use joins more. They are actually more intuitive.
It's not that relational algebra is untintuitive. It's because standard SQL sucks.
Random order is also a valid order.
People who want it, want it with an order.
Look at
https://stackoverflow.com/questions/121387/fetch-the-row-whi...
and
https://stackoverflow.com/questions/3800551/select-first-row...
and
https://stackoverflow.com/questions/8748986/get-records-with...
and their combined thousands of votes and dozens of answers, all full of awkward workaround or ill-performing or specialised-for-one-database-engine code for this common and desirable thing which would be trivial with a couple of boring loops in Python.
I was thinking SQL could do with a keyword for that, like maybe FILTER, that looks like a JOIN but works like EXISTS.
[OUTER|SEMI|ANTI|ANY|ASOF]
WHERE x IN (SELECT ...)
Maybe I am missing something, but the exact meaning to vary a lot across different databases. Some seem to have a standalone "table" data type, while others don't.
How it's implemented is as usual up to the database server implementation. Ones I've used creates a temporary table (like it does in so many other cases), and as such EXISTS is usually faster.
But I wouldn't rely on this when moving to another implementation, and use the query planner to see, just as I'd view the assembly output when moving to a new compiler.
Again, I don't have tons of experience, so concrete (counter) examples are welcome.
In my experience, in the non-ad-hoc use-case, views can often be substituted for the procedural approach, forming the equivalent of a POSIX pipe.
*> A colleague of mine wrote a paper many years ago, pointing out that thinking procedurally is more "natural" for many: https://dl.acm.org/doi/10.1145/319628.319656. But thinking set-at-a-time instead of one-at-a-time is a valuable skill, not that far off from thinking functionally.
Hmm. Given the proliferation of tabular data tools (especially spreadsheets) over the intervening 40 years, I wonder if those results would remain the same today (and whether there would be any difference among Excel power users that use pivot tables, etc.)