I didn’t really see loops being handled here. As far as I understand, the biggest technical difficulty with this kind of probabilistic programming language is handling loops, including infinite loops and almost surely terminating loops.
I did a bit of research earlier in my life[0] to study the handling of loops and without using Monte Carlo simulation. The result was actually workable if incredibly resource intensive to the point of being impractical. If I had chosen to do it again, I might’ve accepted using Monte Carlo simulations while still supporting loops.
I started this about 9 years ago and never finished it. The idea comes from a course in my telecom degree called "Señales Aleatorias y Ruido" (Random Signals and Noise), I spent so many evenings writing probability by hand, and every time I wanted to check a result with a computer it was a ton of boilerplate.
The engine is Rust, the JIT is built on Cranelift, there is also a WASM backend so everything runs in the browser too.
Full disclosure, I could only finish it now because of AI agents. In my experience they are amazing at the runtime and the numerical code, but pretty bad at language design, so I kept that part for myself.
I built a small one of these as a calculator back in the day, it was fun and I always wanted to see someone really run with the idea more, so your project is super cool!
Did you ever look at symbolic or exact operators instead of purely monte-Carlo?
I remember reading the paper for distr, they used a mix of Fourier transforms for convolution and symbolic reduction to build a probabilistic computing library in R. I attempted building a small python library for this, but for my problems the CLT ended up sufficient to approximate the results faster, so I went with that.
You may enjoy reading the paper, it’s not groundbreaking but is a nice presentation of relationships/operations. Maybe it’ll inspire some features for you.
1. It's obvious that you are a big fan of Cranelift. I'd be interested to hear more about your experience in practical terms. For example could you share any insight about use cases where it is best suited, and where it might be better to look elsewhere? Did you hit any pain points? What was its killer feature for NoiseLang?
2. You wrote: "My favorite trick is in the RNG. Generating random numbers is a serial dependency chain, so instead of fighting that, the kernel runs four independent streams at once and lets the out-of-order core overlap them. This trick ended up beating a hand-written SIMD kernel!" What does this mean exactly? you just ran a scalar kernel 4-wide using SIMD instructions? or you interleaved 4 scalar copies of the same algorithm? did you generate the code or just duplicate the streams by hand?
I was exploring ways to speed up this language, the naive implementation is just a interpreter executed in rust, which can just do so much. Once thing i explored was to compile the program graph into WASM, then execute WASM. The idea is the the WASM runtime would JIT program and run faster than any interpreter I could write myself. During this exploration, i found that I could use the JIT optimizer directly and skip the WASM step.
What noiselang does it parse, convert to a execution graph then carefully use the Cranelift api to describe the program, and under the hood, it will find different opetimization and generate byte code that runs directly in the host CPU.
The reason for it to be killer is that it allows NoiseLang to run as native speeds, with very little compiler/optimization work. it's a very simple repo.
For the RNG, this was a discovery myself, when profiling, i found the many benchmarks were limited by the speed of the RNG itself, ie, if i could genenrate random numbers faster, the simulation would be faster. xoshiro's next number is computed from its current state. So to get number N+1 you must have finished number N. It's a chain: A → B → C → D. Your CPU can run maybe 6 integer operations per cycle, but a chain only ever offers it one to run. Five of the six lanes sit empty.
I tried to use SIMD to speed this up, but still hit the limit, even if using SIMD, it still had to wait for the next number, a massive speed up came from realizing that i can keep four independent xoshiro256++ states and emits four samples per loop iteration, i += 4. Since the four state-update chains share no registers, the out-of-order core issues them in the same cycles instead of stalling on one serial chain.
SIMD gives you more work per instruction. But I wasn't short on work, I was waiting. A 2-wide xoshiro still needs state N before it can compute state N+1, so the chain is the same length and I wait at every link, I just get two numbers per link instead of one.
And each link costs more. xoshiro rotates a 64-bit word every round, and NEON has no 64-bit rotate, so that becomes three instructions instead of one. Twice the numbers, three times the wait.
Four streams wins because it leaves the chain alone. It just runs four of them at once, and the CPU was already idle enough to overlap them for free.
Surely SIMD combined with multiple streams would beat both approaches. (This would be separate streams in each SIMD lane and separate streams in different SIMD variables.) There are multiple SIMD execution units, just like the 6 scalar units you mention. The latency of SIMD ops will be similar to scalar, except in cases you mention like shifts.
oh! that is very interesting. I was not aware of I could simulate markov chains with Approximate Bayesian, I have some good reading to do this weekend! indeed, expressions like P(D == 8 | D > 3) are already natively supported: https://noiselang.com/play/#x=conditional_bayes
Fair! My thinking was that PM of a single tone signal (the one i use in the demo is equivalent to FM, but shifted a bit). And implementing real FM for decoding is a lot more noisy, but I will add some callout in the article.
Truth be told, you motivated me to write the exact FM with the differenciation, maybe. Could be interesting to simulate PM vs FM for non single tone signals, to see how FM does even better!
I know haha I wanted to be transparent about this, I have been coding since 9 years old, 32 years old now. I have nothing to prove other than it would have been impossible for me find time to complete this project without help, also a Toy language. Not trying to replace anything people use today :) it's a cute project
Since the advent of AI I have been delighted by the number of old project ideas I've been able to execute. There are just too many ideas to implement them all by yourself - but AIs don't seem to mind in the least.
Claude Max Plan honestly is kind of endless, I have a old game written in MRC Objectice-C that i am porting to something more modern, and still within the limits of the flat subscription... not sure for how long this will last
Interestingly, shading languages started out like this - way before consumer GPUs.
I remember encountering this idea written in a book written by Ed Catmull of Pixar fame (can't find the title sorry, but it was written in the 80s), but generally comes from signal processing as a way of avoiding aliasing artifacts..
The core idea is to make programming, which is a discrete and discontinuous domain, into a well-behaved band limited signal. Otherwise you get aliasing (or jaggies), which can happen even INSIDE a surface, if the shader's like that.
The code idea for this is the step function which is the integral of the dirac delta. step(x) returns 1 for all x >0 and 0 otherwise.
Step is not a well-behaved function in the sense, that it changes infinitely quickly at x=0. But once we know what we want, we can replace it with something like that, that's well behaved.
Consider the example pseudocode
color = x> 5? green:blue;
can be rewritten as
color = blue + step(x-5)(green-blue)
With the two being equivalent.
Now if we put the code into a shader, we get jaggies. So to combat the value changing infinitely fast, we go for a function that's like step, but changes smoothly* from 0 to 1 around x=0. Enter smoothstep:
color = blue + smoothstep(x-(5+EPSILION),(x-EPSILON), x)*(green-blue)
And so we defined a 'transition zone' of +-EPSILON(an arbitrary number). While any smooth function can work, smoothstep is chosen because it has a smooth first and second derivative (meaning even if you want to get the rate of change, something that often pops up in computer graphics, the result will be still well behaved).
Pixar's Renderman shading language (which is remarkably similar to GLSL/HLSL/C), used to do this automatically for you. Essentially it could take arbitrary code peppered with if statements, and turn it into a continuous function.
Which is kinda cool imo.
It's also a cool trick in the age of AI. Since you have a function that's well-behaved, you can do things like gradient descent to train an AI to synthetize a function for you. You can even say, that you don't need exact results, you can accept some error.
In this case your program optimization problem can be reframed from doing idempotent transformations on the list of instructions, to getting a program that generates a target function whose error is no greater than some (mathematical) reference function.
this is also related somewhat to the notion of differentiable programming. RELU is (roughly) the same as x * step(x). In differentiable programming one can replace it with smooth approximations, cf "softplus"
That book also has a chapter on control flow, which is very similar to what you're talking about.
Unrolling an if statement into x = b (result of one branch) + (1-b) (result of the other branch) is also incredibly common in cryptography. If `b` is a "secret" variable, an if statement may leak the value of it via the branch predictor/speculative execution. The way around this is to compute both branches, and then select them with the above arithmetic expression. This mostly works, though compilers are tediously smart, and so one often has to be careful how with how you precisely do it.
Nice! I’ve dabbled with something similar on my own lately (originally wrote/vibed to explain some concepts that came up when discussing D&D) at diceplots.com - different approach, keeping the distributions exactly analytical at every step, never sampling.
Does it still count as a Dirac delta when it’s a discrete distribution? (The distributions in TFA are not continuous - they are things like a roll of 1d6 etc)
Yes, a Dirac delta is just "all the weight on one point", and that works fine on a die.
For the scope of the language it never even comes up, because Noise is a simulator, it does not evaluate densities, it draws samples.
The point is that every value goes through the same operators. Add them, compare them, pass them to a function, put one in the condition of an if. You can even use a random variable to define another random variable:
bias ~ unif(0, 1)
flips ~[10] bernoulli(bias) // bernoulli just took a distribution where a number normally goes.
But you right, dirac only applies to continuous functions, in Noise is only refers to the dirac measure. I found this article a fun/nerd to make my point that everything "acts" as a distribution from the DX perspective, but under the hood 5 is just 5.
And a constant collapses back to a plain integer in the graph anyway, so 5 costs nothing.
I like the “everything is a distribution” approach a lot[1]. Looks great. I’m looking forward to actually having a chance to mess around with it.
[1] And feels philosophically like the unification in the underlying maths between discrete and continuous probability that you get when you apply measure theory
Stan and PyMC beat Noise at the thing they’re built for, fitting a posterior to lots of continuous data with their HMC/NUTS samplers, and NumPy beats it at raw array crunching. Conditioning in Noise is rejection-based, so it works great for a handful of discrete observations but becomes useless for ten thousand continuous measurements, and there is no stateful simulation yet (no Markov chains yet). Where Noise wins when you have a probability question and you wanna know the answer without much hassle.
So use Noise for the whiteboard stage of a problem, when you want to run the math you just wrote, and move to Stan or PyMC when you need a real posterior, or to NumPy and JAX when you need to go to production.
People do this kind of computation using numpy all the time. If you weren’t developing a new language with a new syntax, there really isn’t much of a library to write. At that point it’s just using numpy.
the language can do a lot of very expressive things, every language feature works in your favor too, but agree with you, i would not use my own language for anything production.
like this:
D ~ unif_int(1, 6);
Print("P(rolled a 6 | rolled > 3) =", P(D == 6 | D > 3));
or:
loss ~ unif(0, 1000);
claim = if loss > 200 { loss - 200 } else { 0 };
p = P(claim > 0);
Print("P(insurer pays a claim) =", p)
notice that "claim" is also a random variable! result of a if expression
N = 10**6
D = np.random.randint(1, 6, N)
print("P(rolled a 6 | rolled > 3) =", ((D == 6) | (D > 3)).mean())
loss = np.random.uniform(0, 1000, N)
claim = np.where(loss > 200, loss - 200, 0)
p = (claim > 0).mean()
print("P(insurer pays a claim) =", p)
It’s concise enough that people generally wouldn’t bother writing a library. Unless they really want their custom syntax, then perhaps they write a parser.
that has some technical limitations. For example, their impl can compile to wasm, which makes giving an online interpreter simpler/lighter weight than relying on running python in the browser.
My system is blocking that site as it is on the HaGeZi blocklist. I don't have any further information, and I'm not expressing an opinion on the site. An alternative might be https://noiselang.com, which is not on the blocklist.
mmmh i can't see the domain blocked in the list, it's my personal blog, i don't even have tracking other than server-side stats. could it be because using netlify dns?
Firstly, I'm not intending any slight on you personally! In fact this might be more of an issue for you interacting with the site than for people just reading an article.
I didn’t really see loops being handled here. As far as I understand, the biggest technical difficulty with this kind of probabilistic programming language is handling loops, including infinite loops and almost surely terminating loops.
I did a bit of research earlier in my life[0] to study the handling of loops and without using Monte Carlo simulation. The result was actually workable if incredibly resource intensive to the point of being impractical. If I had chosen to do it again, I might’ve accepted using Monte Carlo simulations while still supporting loops.
[0]: https://github.com/kccqzy/probabilistic-program-inference/bl... Shameless self promotion I know! I put quite a bit of effort into that README and the code.
I started this about 9 years ago and never finished it. The idea comes from a course in my telecom degree called "Señales Aleatorias y Ruido" (Random Signals and Noise), I spent so many evenings writing probability by hand, and every time I wanted to check a result with a computer it was a ton of boilerplate.
The engine is Rust, the JIT is built on Cranelift, there is also a WASM backend so everything runs in the browser too.
Full disclosure, I could only finish it now because of AI agents. In my experience they are amazing at the runtime and the numerical code, but pretty bad at language design, so I kept that part for myself.
It's a toy language. Ask me anything!
I built a small one of these as a calculator back in the day, it was fun and I always wanted to see someone really run with the idea more, so your project is super cool!
Did you ever look at symbolic or exact operators instead of purely monte-Carlo?
I remember reading the paper for distr, they used a mix of Fourier transforms for convolution and symbolic reduction to build a probabilistic computing library in R. I attempted building a small python library for this, but for my problems the CLT ended up sufficient to approximate the results faster, so I went with that.
You may enjoy reading the paper, it’s not groundbreaking but is a nice presentation of relationships/operations. Maybe it’ll inspire some features for you.
https://arxiv.org/pdf/1006.0764
Very cool. I have a couple of questions:
1. It's obvious that you are a big fan of Cranelift. I'd be interested to hear more about your experience in practical terms. For example could you share any insight about use cases where it is best suited, and where it might be better to look elsewhere? Did you hit any pain points? What was its killer feature for NoiseLang?
2. You wrote: "My favorite trick is in the RNG. Generating random numbers is a serial dependency chain, so instead of fighting that, the kernel runs four independent streams at once and lets the out-of-order core overlap them. This trick ended up beating a hand-written SIMD kernel!" What does this mean exactly? you just ran a scalar kernel 4-wide using SIMD instructions? or you interleaved 4 scalar copies of the same algorithm? did you generate the code or just duplicate the streams by hand?
I was exploring ways to speed up this language, the naive implementation is just a interpreter executed in rust, which can just do so much. Once thing i explored was to compile the program graph into WASM, then execute WASM. The idea is the the WASM runtime would JIT program and run faster than any interpreter I could write myself. During this exploration, i found that I could use the JIT optimizer directly and skip the WASM step.
What noiselang does it parse, convert to a execution graph then carefully use the Cranelift api to describe the program, and under the hood, it will find different opetimization and generate byte code that runs directly in the host CPU.
The reason for it to be killer is that it allows NoiseLang to run as native speeds, with very little compiler/optimization work. it's a very simple repo.
For the RNG, this was a discovery myself, when profiling, i found the many benchmarks were limited by the speed of the RNG itself, ie, if i could genenrate random numbers faster, the simulation would be faster. xoshiro's next number is computed from its current state. So to get number N+1 you must have finished number N. It's a chain: A → B → C → D. Your CPU can run maybe 6 integer operations per cycle, but a chain only ever offers it one to run. Five of the six lanes sit empty.
I tried to use SIMD to speed this up, but still hit the limit, even if using SIMD, it still had to wait for the next number, a massive speed up came from realizing that i can keep four independent xoshiro256++ states and emits four samples per loop iteration, i += 4. Since the four state-update chains share no registers, the out-of-order core issues them in the same cycles instead of stalling on one serial chain.
SIMD gives you more work per instruction. But I wasn't short on work, I was waiting. A 2-wide xoshiro still needs state N before it can compute state N+1, so the chain is the same length and I wait at every link, I just get two numbers per link instead of one.
And each link costs more. xoshiro rotates a 64-bit word every round, and NEON has no 64-bit rotate, so that becomes three instructions instead of one. Twice the numbers, three times the wait.
Four streams wins because it leaves the chain alone. It just runs four of them at once, and the CPU was already idle enough to overlap them for free.
Surely SIMD combined with multiple streams would beat both approaches. (This would be separate streams in each SIMD lane and separate streams in different SIMD variables.) There are multiple SIMD execution units, just like the 6 scalar units you mention. The latency of SIMD ops will be similar to scalar, except in cases you mention like shifts.
Definitely going to play around with this, thanks for posting.
I know MCMC isn’t your goal, but seems like this could be used for ABC-MCMC (as is?)
Would also be nice to have an option to plot using a KDE vs histograms.
(Also your FM example seems to be technically PM)
oh! that is very interesting. I was not aware of I could simulate markov chains with Approximate Bayesian, I have some good reading to do this weekend! indeed, expressions like P(D == 8 | D > 3) are already natively supported: https://noiselang.com/play/#x=conditional_bayes
Fair! My thinking was that PM of a single tone signal (the one i use in the demo is equivalent to FM, but shifted a bit). And implementing real FM for decoding is a lot more noisy, but I will add some callout in the article.
Truth be told, you motivated me to write the exact FM with the differenciation, maybe. Could be interesting to simulate PM vs FM for non single tone signals, to see how FM does even better!
Be warned - by using AI like this you've made yourself a lightning rod for the people who really really really dislike AI.
I know haha I wanted to be transparent about this, I have been coding since 9 years old, 32 years old now. I have nothing to prove other than it would have been impossible for me find time to complete this project without help, also a Toy language. Not trying to replace anything people use today :) it's a cute project
Yeah. :)
Since the advent of AI I have been delighted by the number of old project ideas I've been able to execute. There are just too many ideas to implement them all by yourself - but AIs don't seem to mind in the least.
It's a brave new world.
how much did Ai cost for you?
Claude Max Plan honestly is kind of endless, I have a old game written in MRC Objectice-C that i am porting to something more modern, and still within the limits of the flat subscription... not sure for how long this will last
so 200$
Interestingly, shading languages started out like this - way before consumer GPUs.
I remember encountering this idea written in a book written by Ed Catmull of Pixar fame (can't find the title sorry, but it was written in the 80s), but generally comes from signal processing as a way of avoiding aliasing artifacts..
The core idea is to make programming, which is a discrete and discontinuous domain, into a well-behaved band limited signal. Otherwise you get aliasing (or jaggies), which can happen even INSIDE a surface, if the shader's like that.
The code idea for this is the step function which is the integral of the dirac delta. step(x) returns 1 for all x >0 and 0 otherwise. Step is not a well-behaved function in the sense, that it changes infinitely quickly at x=0. But once we know what we want, we can replace it with something like that, that's well behaved.
Consider the example pseudocode
can be rewritten as color = blue + step(x-5)(green-blue)With the two being equivalent.
Now if we put the code into a shader, we get jaggies. So to combat the value changing infinitely fast, we go for a function that's like step, but changes smoothly* from 0 to 1 around x=0. Enter smoothstep: color = blue + smoothstep(x-(5+EPSILION),(x-EPSILON), x)*(green-blue)
And so we defined a 'transition zone' of +-EPSILON(an arbitrary number). While any smooth function can work, smoothstep is chosen because it has a smooth first and second derivative (meaning even if you want to get the rate of change, something that often pops up in computer graphics, the result will be still well behaved).
Pixar's Renderman shading language (which is remarkably similar to GLSL/HLSL/C), used to do this automatically for you. Essentially it could take arbitrary code peppered with if statements, and turn it into a continuous function.
Which is kinda cool imo.
It's also a cool trick in the age of AI. Since you have a function that's well-behaved, you can do things like gradient descent to train an AI to synthetize a function for you. You can even say, that you don't need exact results, you can accept some error.
In this case your program optimization problem can be reframed from doing idempotent transformations on the list of instructions, to getting a program that generates a target function whose error is no greater than some (mathematical) reference function.
note that similar concepts appear in mathematics. Generally the term for it is a "mollified" function.
applied to the step function, you would get a smooth cutoff function
https://en.wikipedia.org/wiki/Mollifier#Smooth_cutoff_functi...
this is also related somewhat to the notion of differentiable programming. RELU is (roughly) the same as x * step(x). In differentiable programming one can replace it with smooth approximations, cf "softplus"
https://arxiv.org/pdf/2403.14606
That book also has a chapter on control flow, which is very similar to what you're talking about.
Unrolling an if statement into x = b (result of one branch) + (1-b) (result of the other branch) is also incredibly common in cryptography. If `b` is a "secret" variable, an if statement may leak the value of it via the branch predictor/speculative execution. The way around this is to compute both branches, and then select them with the above arithmetic expression. This mostly works, though compilers are tediously smart, and so one often has to be careful how with how you precisely do it.
I don't know JAX, but can this trick be applied as a higher-level function to autodiff-capable languages like JAX?
Reminds me of Haskell’s monad-bayes: https://monad-bayes.netlify.app/
oh! that's awesome, i had no idea haskell could express this things
Nice! I’ve dabbled with something similar on my own lately (originally wrote/vibed to explain some concepts that came up when discussing D&D) at diceplots.com - different approach, keeping the distributions exactly analytical at every step, never sampling.
Does it still count as a Dirac delta when it’s a discrete distribution? (The distributions in TFA are not continuous - they are things like a roll of 1d6 etc)
Yes, a Dirac delta is just "all the weight on one point", and that works fine on a die.
For the scope of the language it never even comes up, because Noise is a simulator, it does not evaluate densities, it draws samples.
The point is that every value goes through the same operators. Add them, compare them, pass them to a function, put one in the condition of an if. You can even use a random variable to define another random variable:
bias ~ unif(0, 1) flips ~[10] bernoulli(bias) // bernoulli just took a distribution where a number normally goes.
and in if-stataments:
DistributionC = if DistributionA < DistributionB { 0 } else { 1 }
But you right, dirac only applies to continuous functions, in Noise is only refers to the dirac measure. I found this article a fun/nerd to make my point that everything "acts" as a distribution from the DX perspective, but under the hood 5 is just 5.
And a constant collapses back to a plain integer in the graph anyway, so 5 costs nothing.
I like the “everything is a distribution” approach a lot[1]. Looks great. I’m looking forward to actually having a chance to mess around with it.
[1] And feels philosophically like the unification in the underlying maths between discrete and continuous probability that you get when you apply measure theory
It might be worth looking into probabilistic programming languages. I'm out of date, but I remember webppl, stan, anglican, pymc (a python library).
Seems worth an investigation and maybe mention on the article.
if you read the website, the author explicitly mentions stan in the comparison at the end ;^)
It says
Stan and PyMC beat Noise at the thing they’re built for, fitting a posterior to lots of continuous data with their HMC/NUTS samplers, and NumPy beats it at raw array crunching. Conditioning in Noise is rejection-based, so it works great for a handful of discrete observations but becomes useless for ten thousand continuous measurements, and there is no stateful simulation yet (no Markov chains yet). Where Noise wins when you have a probability question and you wanna know the answer without much hassle.
So use Noise for the whiteboard stage of a problem, when you want to run the math you just wrote, and move to Stan or PyMC when you need a real posterior, or to NumPy and JAX when you need to go to production.
This reminds me of https://mc-stan.org
You should read the part of the article where the author compares his work to Stan.
I'd have just written this as a Python library that lazily evaluates expression via numpy personally. The API is useful, language is not
People do this kind of computation using numpy all the time. If you weren’t developing a new language with a new syntax, there really isn’t much of a library to write. At that point it’s just using numpy.
the language can do a lot of very expressive things, every language feature works in your favor too, but agree with you, i would not use my own language for anything production.
like this:
D ~ unif_int(1, 6); Print("P(rolled a 6 | rolled > 3) =", P(D == 6 | D > 3));
or:
loss ~ unif(0, 1000); claim = if loss > 200 { loss - 200 } else { 0 }; p = P(claim > 0); Print("P(insurer pays a claim) =", p)
notice that "claim" is also a random variable! result of a if expression
Yeah sure this is how numpy people would do it:
It’s concise enough that people generally wouldn’t bother writing a library. Unless they really want their custom syntax, then perhaps they write a parser.that has some technical limitations. For example, their impl can compile to wasm, which makes giving an online interpreter simpler/lighter weight than relying on running python in the browser.
WOWW!!!!!
My system is blocking that site as it is on the HaGeZi blocklist. I don't have any further information, and I'm not expressing an opinion on the site. An alternative might be https://noiselang.com, which is not on the blocklist.
mmmh i can't see the domain blocked in the list, it's my personal blog, i don't even have tracking other than server-side stats. could it be because using netlify dns?
Firstly, I'm not intending any slight on you personally! In fact this might be more of an issue for you interacting with the site than for people just reading an article.
There are multiple versions of the list. The authoritative site appears to be https://github.com/hagezi/dns-blocklists, and making a fairly random choice, I used the "medium" version of the "Threat intelligence feed", and specifically the one marked "Link" for AdBlock. That took me to https://cdn.jsdelivr.net/gh/hagezi/dns-blocklists@latest/adb..., and manualmeida.dev does appear in that list.
The software I'm using is Little Snitch on a Mac, but since the entry is in the list, that's not the problem.
Sent issue: https://github.com/hagezi/dns-blocklists/issues/10759
This is great! thank you for the digging, gonna figuring out how it end up there