8 ms·

When Greedy Algorithms Can Be Faster [C++]

63 points by def-pri-pub 2y ago

cozzyd 2y ago
I suspect if you started using SIMD instructions, the analytical case would get better again (since it's branchless).
- nxobject 2y ago
  Apropos of SIMD – I'm also surprised that the inner loop of the rejection-based algorithm optimized to MMX, but not the analytic algorithm! I would like to think the rejection algorithm after -O3 is benefiting from branch prediction and all sorts of modern speculation optimizations. But I imagine the real test of that would be running these benchmarks would be running these benchmarks on a 5-10ish year old uarch.
  - Sesse__ 2y ago
    Ten years ago, you already have Skylake with pretty good indirect branch predictors...
    - Veliladon 2y ago
      It doesn't matter because you can't predict random data which is what the example is using. Branch predictors work on the premise that 95% of the time the branch resolves the same way which is damn common in most code. The closer your branches look to a normal distribution vs bimodal, the worse the branch predictor is going to perform. The rejector is going to be faster because more iterations can be kept in the reorder buffer at once. The CPU is going to have to keep all the instructions post-branch in the reorder buffer until the branch is retired. If it's waiting on an instruction with an especially long latency like a square root it's probably speculatively already finished the work anyway, rejected or not. If the branch rejects after the work is done the reorder buffer can retire each instruction of the random generation algorithm as it comes and then wait on however many branches at the end which can also be run in parallel because those branches aren't dependent on each other. All those branches which are doing a square root will also be pipelined properly instead of putting bubbles everywhere.
      - Sesse__ 2y ago
        My point was that there hasn't been a revolution in branch prediction (or really any sort of speculation) the last 5–10 years. Things have become incrementally better as always, and the M1 presumably has some sort of pointer value speculation, but by and large, things are generally similar. So running on a 5–10 year old microarchitecture won't change the qualitative results much. If you talked about running on a 30-year-old architecture, then sure, that would tell you something about the effect of modern speculation on the code.
- shihab 2y ago
  On GPU, the rejection sampling approach wouldn't come close to analytical one.
- camel-cdr 2y ago
  With modern SIMD the rejection method is likely faster, because you can do it branchless with a compressed store.
nxobject 2y ago
I don't know how numerics in hardware works, but would the use of functions like sin, cos, sqrt incur a penalty as well, even if only a slight one? It's really fascinating to think about how all of this would work.
- hinkley 2y ago
  Very early games had lookup tables for trig functions. The cpu instructions were too slow or missing. The tables were either generated at run time or statically defined in the code. I think that’s one of those things Jai and Zig agree on - compile time functions have a place in preventing magic numbers that cannot be debugged.
- shihab 2y ago
  yeah, that's very likely the explanation. All these functions are pretty high latency instructions, vs rejection sampling which only involves a multiplication. On Nvidia GPUs, mul has latency of 1-4 cycles while others are 16-32.
- akoboldfrying 2y ago
  Yes, and at least for sqrt(), internally it's likely implemented as a heuristic guess followed by a fixed number of iterations of Newton's Method. (In software, you'd normally iterate Newton's Method until the change in the result is less than some threshold; in hardware, I'm guessing that it might be simpler to figure out the maximum number of iterations that would ever be needed for any input, and always run that many, but I don't know.)
  - Const-me 2y ago
    > at least for sqrt(), internally it's likely implemented as a heuristic guess Square roots are implemented in hardware: https://www.felixcloutier.com/x86/sqrtsd https://www.felixcloutier.com/x86/sqrtsd > In software, you'd normally iterate Newton's Method Software normally computes trigonometric functions (and other complicated ones like exponents and std::erf) with a high-degree polynomial approximation.
Negitivefrags 2y ago
This isn't really what the article is about, but I don't think that the term "Greedy algorithm" means what the author thinks. Greedy algorithms are about making locally optimal choices. They are not "brute force" algorithms or inefficient ones. In fact, greedy algorithms are almost always faster. They are faster because they consider only local information instead of the entire data set. In exchange for that, a greedy algorithm may produce non-optimal answers.
- adamgordonbell 2y ago
  Yeah, greedy is not what he thinks it is. I think of change making algorithm when working retail. Greedy is you always use the largest coin you can, and then the next largest and so on. It works with sensible coin denominations but there are sets of coins where the greedy algorithm is not optimal.
  - o11c 2y ago
    Relevant link: https://en.wikipedia.org/wiki/Change-making_problem https://en.wikipedia.org/wiki/Change-making_problem Looking at cases where greedy isn't optimal, I see two patterns. Either: * there are two coins close in value (ratio of less than 2), e.g. both 20¢ and 25¢; or * there is a "missing" coin at the GCD of two larger coins. I'm pretty sure you can precompute extra "virtual coins" (e.g. 40¢) to make greedy optimal again, but you have to place restrictions on how many of them you're allowed to use.
Labo333 2y ago
To accelerate the analytical solution, there are more obvious techniques, like replacing trigonometric opetations with series. This is a 10x improvement in https://ashvardanian.com/posts/google-benchmark/#slow-trigonometry-in-libc-and-stl-and-fast-math https://ashvardanian.com/posts/google-benchmark/#slow-trigon...
emtel 2y ago
I’m surprised they didn’t mention the reason that the rejection sampling method is surprisingly fast: the probability of needing to resample N times decreases exponentially with N. So even for the 3D case, where 50% of your samples get rejected, the EV for number of samples required is about 2. This is also a good case study on the difference between throughput sensitive vs latency sensitive applications. The rejection sampling is fast on average, but the analytical solution likely has a much tighter upper bound on how long it can take. In a throughput application you care about EV. But if you need to minimize latency, then you care about the upper bound.
chmod775 2y ago
Isn't there a likely faster analytical solution? At least I think sqrt, sin, and cos are generally considered to be slow. Maybe this: 1. Generate random X [-1 to 1] 2. Based on X you can calculate a range of legal Y that will still fall within the circle. Try to use a lookup table+interpolation with desired resolution to make this really fast if there's not some sufficiently fast hardware instructions we can abuse. 3. Generate random Y in that range.
- taco9999 2y ago
  A potential issue with that is needing to adjust the distribution of either X or Y, since if X and Y remain uniform, then there will be a lot of points packed into the sides of the circle with higher density than the center column.
- Sesse__ 2y ago
  > At least I think sqrt, sin, and cos are generally considered to be slow. […] Based on X you can calculate a range of legal Y that will still fall within the circle. Yup, that calculation is called sqrt(1 - x²). > Try to use a lookup table+interpolation with desired resolution to make this really fast On modern hardware (like, most architectures from 1995-ish onwards), LUTs are mostly slower than polynomial-based solutions for the same level of accuracy. > 1. Generate random X [-1 to 1] As others have pointed out, this is already wrong unless you intend to reject some of the solutions later.
psanchez 2y ago
I think I would call it naive algorithm rather than greedy. It looked like an interesting problem so I spent some time this morning exploring if there would be any performance improvement by pregenerating an array of X items (where X around 1M to 16M items) and then randomly returning one of them at a time. I explored the project and copied the functions to be as faithful as the original implementation as possible. Generating 10M unit sphere (best of 3 runs, g++ 13, linux, Intel i7-8565U, one core for tests): - naive/rejection: ~574ms - analytical: ~1122ms - pregen 1M elements: ~96ms That's almost 6x faster than rejection method. Setup of the 1M elements is done once and does not count on the metrics. Using double type, using float yields around 4x improvements. After looking at those results I decided to try on the project itself, so I downloaded, compiled and applied similar optimizations in the project, only updating circle and sphere random generators (with 16M unit vectors that are only created once on app lifetime) but got almost no noticeable benefits (marginal at most). Hard to tell because of the random nature of the raytracing implementation. On the bright side the image quality was on par. Honestly I was afraid this method would generate poor visuals. Just for the record, I'm talking about something as simple as: std::vector<Vec3> g_unitSpherePregen; uint32_t g_unitSpherePregenIndex = 0; void setupUnitSpherePregen(uint32_t nElements) { g_unitSpherePregen.resize(nElements); for (auto i = 0; i < nElements; i++) { g_unitSpherePregen[i] = unitSphereNaive(); // call the original naive or analytical method } } Vec3 unitSpherePregen() { g_unitSpherePregenIndex = (g_unitSpherePregenIndex + 1) % g_unitSpherePregen.size(); return g_unitSpherePregen[g_unitSpherePregenIndex]; } I tried as well using a psrng (std::mt19937 and xorshf96) in unitSpherePregen instead of the incremented variable, but increment was faster and yielded good visual results. Next step would be profiling, but I don't think I will invest more time on this. Edit: fix formatting
- camel-cdr 2y ago
  I also came up with an alternative implementation: https://gist.github.com/camel-cdr/d16fd2be1fd7b71622649e6bc7fd224a https://gist.github.com/camel-cdr/d16fd2be1fd7b71622649e6bc7... The idea is based on the Ziggurat Method. You overlap the circle with n boxes that each encapsulate the same amount of area of the underlying circle, select a random box, and then do rejection. With 128 boxes, this reduces the average number of additional iterations from 27% to 0.7%, which should massively reduce the number of branch miss predictions. It ended up about 2x faster the simple rejection method. I haven't applied this to spheres yet, but that should also be possible.

63 points by def-pri-pub 2y ago

cozzyd 2y ago

I suspect if you started using SIMD instructions, the analytical case would get better again (since it's branchless).

nxobject 2y ago
Apropos of SIMD – I'm also surprised that the inner loop of the rejection-based algorithm optimized to MMX, but not the analytic algorithm! I would like to think the rejection algorithm after -O3 is benefiting from branch prediction and all sorts of modern speculation optimizations. But I imagine the real test of that would be running these benchmarks would be running these benchmarks on a 5-10ish year old uarch.
- Sesse__ 2y ago
  Ten years ago, you already have Skylake with pretty good indirect branch predictors...
  - Veliladon 2y ago
    It doesn't matter because you can't predict random data which is what the example is using. Branch predictors work on the premise that 95% of the time the branch resolves the same way which is damn common in most code. The closer your branches look to a normal distribution vs bimodal, the worse the branch predictor is going to perform. The rejector is going to be faster because more iterations can be kept in the reorder buffer at once. The CPU is going to have to keep all the instructions post-branch in the reorder buffer until the branch is retired. If it's waiting on an instruction with an especially long latency like a square root it's probably speculatively already finished the work anyway, rejected or not. If the branch rejects after the work is done the reorder buffer can retire each instruction of the random generation algorithm as it comes and then wait on however many branches at the end which can also be run in parallel because those branches aren't dependent on each other. All those branches which are doing a square root will also be pipelined properly instead of putting bubbles everywhere.
    - Sesse__ 2y ago
      My point was that there hasn't been a revolution in branch prediction (or really any sort of speculation) the last 5–10 years. Things have become incrementally better as always, and the M1 presumably has some sort of pointer value speculation, but by and large, things are generally similar. So running on a 5–10 year old microarchitecture won't change the qualitative results much. If you talked about running on a 30-year-old architecture, then sure, that would tell you something about the effect of modern speculation on the code.
shihab 2y ago
On GPU, the rejection sampling approach wouldn't come close to analytical one.
camel-cdr 2y ago
With modern SIMD the rejection method is likely faster, because you can do it branchless with a compressed store.

nxobject 2y ago

I don't know how numerics in hardware works, but would the use of functions like sin, cos, sqrt incur a penalty as well, even if only a slight one? It's really fascinating to think about how all of this would work.

hinkley 2y ago
Very early games had lookup tables for trig functions. The cpu instructions were too slow or missing. The tables were either generated at run time or statically defined in the code. I think that’s one of those things Jai and Zig agree on - compile time functions have a place in preventing magic numbers that cannot be debugged.
shihab 2y ago
yeah, that's very likely the explanation. All these functions are pretty high latency instructions, vs rejection sampling which only involves a multiplication. On Nvidia GPUs, mul has latency of 1-4 cycles while others are 16-32.
akoboldfrying 2y ago
Yes, and at least for sqrt(), internally it's likely implemented as a heuristic guess followed by a fixed number of iterations of Newton's Method. (In software, you'd normally iterate Newton's Method until the change in the result is less than some threshold; in hardware, I'm guessing that it might be simpler to figure out the maximum number of iterations that would ever be needed for any input, and always run that many, but I don't know.)
- Const-me 2y ago
  > at least for sqrt(), internally it's likely implemented as a heuristic guess Square roots are implemented in hardware: https://www.felixcloutier.com/x86/sqrtsd https://www.felixcloutier.com/x86/sqrtsd > In software, you'd normally iterate Newton's Method Software normally computes trigonometric functions (and other complicated ones like exponents and std::erf) with a high-degree polynomial approximation.

Negitivefrags 2y ago

This isn't really what the article is about, but I don't think that the term "Greedy algorithm" means what the author thinks. Greedy algorithms are about making locally optimal choices. They are not "brute force" algorithms or inefficient ones. In fact, greedy algorithms are almost always faster. They are faster because they consider only local information instead of the entire data set. In exchange for that, a greedy algorithm may produce non-optimal answers.

adamgordonbell 2y ago
Yeah, greedy is not what he thinks it is. I think of change making algorithm when working retail. Greedy is you always use the largest coin you can, and then the next largest and so on. It works with sensible coin denominations but there are sets of coins where the greedy algorithm is not optimal.
- o11c 2y ago
  Relevant link: https://en.wikipedia.org/wiki/Change-making_problem https://en.wikipedia.org/wiki/Change-making_problem Looking at cases where greedy isn't optimal, I see two patterns. Either: * there are two coins close in value (ratio of less than 2), e.g. both 20¢ and 25¢; or * there is a "missing" coin at the GCD of two larger coins. I'm pretty sure you can precompute extra "virtual coins" (e.g. 40¢) to make greedy optimal again, but you have to place restrictions on how many of them you're allowed to use.

Labo333 2y ago

To accelerate the analytical solution, there are more obvious techniques, like replacing trigonometric opetations with series. This is a 10x improvement in https://ashvardanian.com/posts/google-benchmark/#slow-trigonometry-in-libc-and-stl-and-fast-math https://ashvardanian.com/posts/google-benchmark/#slow-trigon...

emtel 2y ago

I’m surprised they didn’t mention the reason that the rejection sampling method is surprisingly fast: the probability of needing to resample N times decreases exponentially with N. So even for the 3D case, where 50% of your samples get rejected, the EV for number of samples required is about 2. This is also a good case study on the difference between throughput sensitive vs latency sensitive applications. The rejection sampling is fast on average, but the analytical solution likely has a much tighter upper bound on how long it can take. In a throughput application you care about EV. But if you need to minimize latency, then you care about the upper bound.

chmod775 2y ago

Isn't there a likely faster analytical solution? At least I think sqrt, sin, and cos are generally considered to be slow. Maybe this: 1. Generate random X [-1 to 1] 2. Based on X you can calculate a range of legal Y that will still fall within the circle. Try to use a lookup table+interpolation with desired resolution to make this really fast if there's not some sufficiently fast hardware instructions we can abuse. 3. Generate random Y in that range.

taco9999 2y ago
A potential issue with that is needing to adjust the distribution of either X or Y, since if X and Y remain uniform, then there will be a lot of points packed into the sides of the circle with higher density than the center column.
Sesse__ 2y ago
> At least I think sqrt, sin, and cos are generally considered to be slow. […] Based on X you can calculate a range of legal Y that will still fall within the circle. Yup, that calculation is called sqrt(1 - x²). > Try to use a lookup table+interpolation with desired resolution to make this really fast On modern hardware (like, most architectures from 1995-ish onwards), LUTs are mostly slower than polynomial-based solutions for the same level of accuracy. > 1. Generate random X [-1 to 1] As others have pointed out, this is already wrong unless you intend to reject some of the solutions later.

psanchez 2y ago

I think I would call it naive algorithm rather than greedy. It looked like an interesting problem so I spent some time this morning exploring if there would be any performance improvement by pregenerating an array of X items (where X around 1M to 16M items) and then randomly returning one of them at a time. I explored the project and copied the functions to be as faithful as the original implementation as possible. Generating 10M unit sphere (best of 3 runs, g++ 13, linux, Intel i7-8565U, one core for tests): - naive/rejection: ~574ms - analytical: ~1122ms - pregen 1M elements: ~96ms That's almost 6x faster than rejection method. Setup of the 1M elements is done once and does not count on the metrics. Using double type, using float yields around 4x improvements. After looking at those results I decided to try on the project itself, so I downloaded, compiled and applied similar optimizations in the project, only updating circle and sphere random generators (with 16M unit vectors that are only created once on app lifetime) but got almost no noticeable benefits (marginal at most). Hard to tell because of the random nature of the raytracing implementation. On the bright side the image quality was on par. Honestly I was afraid this method would generate poor visuals. Just for the record, I'm talking about something as simple as: std::vector<Vec3> g_unitSpherePregen; uint32_t g_unitSpherePregenIndex = 0; void setupUnitSpherePregen(uint32_t nElements) { g_unitSpherePregen.resize(nElements); for (auto i = 0; i < nElements; i++) { g_unitSpherePregen[i] = unitSphereNaive(); // call the original naive or analytical method } } Vec3 unitSpherePregen() { g_unitSpherePregenIndex = (g_unitSpherePregenIndex + 1) % g_unitSpherePregen.size(); return g_unitSpherePregen[g_unitSpherePregenIndex]; } I tried as well using a psrng (std::mt19937 and xorshf96) in unitSpherePregen instead of the incremented variable, but increment was faster and yielded good visual results. Next step would be profiling, but I don't think I will invest more time on this. Edit: fix formatting

camel-cdr 2y ago
I also came up with an alternative implementation: https://gist.github.com/camel-cdr/d16fd2be1fd7b71622649e6bc7fd224a https://gist.github.com/camel-cdr/d16fd2be1fd7b71622649e6bc7... The idea is based on the Ziggurat Method. You overlap the circle with n boxes that each encapsulate the same amount of area of the underlying circle, select a random box, and then do rejection. With 128 boxes, this reduces the average number of additional iterations from 27% to 0.7%, which should massively reduce the number of branch miss predictions. It ended up about 2x faster the simple rejection method. I haven't applied this to spheres yet, but that should also be possible.