Sir, this is a Wendy's ;P

It might be worth posting this on https://math.stackexchange.com/

Excuse me sir, to me all topics here feel like a wendy's so I thought this was the right place. Thanks for the suggestion, I’ll post it there too!

It is definitely the right place!

> Optimizing and doing it with matrices brought it down to 20 minutes

What was the final run time?

11 seconds in colab implementing the basic transformation with reshape, roll, and flatten

4 seconds in google apps script by handling indices in loops instead of rearranging elements

Are you computing a single expectation value, or are you searching for the bet of maximum expected value among all possible bets ?

How long does it take you to compute a single expectation value ?

I have written some single threaded non-simd c++ code for this : https://gist.github.com/unrealwill/d1bc68d1f5c7ee6fe72b76dc5...

When there are 14 matches, there are 3^14 = 4782969 different possible results and the same number of different possible bets.

Summing naively a single expectation takes 0.1s Summing on the hamming neighborhood to compute a single expectation takes 7.6e-5s,

Computing all expectation values and computing the maximum takes 153s (with g++9 but 390s with g++-10 (I don't know why)).

I am not sure whether or not I used the same trick as you. Here are the tricks I use : when probabilities sum to one, the number of iteration to compute a single expectation is not dependent of the number of possible results for a match, because the sum can be factored by grouping all negative results for a single event together by applying the weight as 1-proba.

Some tricks used are representing the result for the 14 matches as a single int32_t and working in base 3. It can probably be even faster, if you work in base 4 instead and replace integer division by bit shifts.

Iterating and in particular iterating combination is excruciatingly slow in python even when using itertools. So do yourself a favor, and write c++ code for these kind of things so that you don't have to have memory allocations inside loops.

Following your suggestion, I’ve implemented the PSA algorithm (the optimized one) in C++. It might be tricky to understand without an explanation, especially since (1) it’s generalized for any dimensionality, and (2) I manually implemented the basic transformation ("oscillating shift") instead of using NumPy’s reshape+roll+flatten, which I believe makes it easier to grasp. But you can verify that it returns the same results as the naive algorithm. For the quiniela case, it shouldn’t take more than 2-3 seconds to compute the sum of neighbors for all points. https://github.com/petopello/PSA/blob/main/PSA.cpp

The key point is that your code takes as input ‘probas’, a 3x14 table with probabilities of independent events. From this, you assume the structure of the space by multiplying result of 14 events. Here, you can apply the trick of grouping negative probabilities. In fact, some of these calculations can be done only once, as many points share the same multiplication path. This is what I implemented in my SHN_BOXES function.

The issue is that our problem isn’t exactly like that. I haven’t gone into too much detail, but you should be able to create a function that takes a pre-built space as input and sums over it, even if you don’t know how the space was built or whether it was constructed with random values.

For details on why this is the case:

Step 1.1: First, calculate, based on a table 3x14 like ‘probas’, which represents how many people have bet on event j for match i, the number of winners in each category if a certain prediction occurs (I assume homogeneity here). This is also a Hamming neighborhood sum, but you can use the boxes algorithm (~1 sec).

Step 1.2: The prizes depend (inversely) on the number of winners. So, apply a formula like this to the previous space: bet_price * coefficients / (winners + bet_price / revenue)

Where 'winners' are the values calculated in 1.1, and the rest are inputs: REVENUE = 1000000.0 PRICE = 0.75 COEFFICIENTS = [0.16, 0.075, 0.075, 0.075, 0.09] //percentaje of revenue correspondy to each category follow game rules

Step 1.3: Now we have the prizes for each prediction. To calculate the value of betting on a specific prediction, we need to do a sum product of the prizes corresponding to that bet (distance <= 4), each multiplied by its probability. So, we multiply the results from the previous formula by the probability of occurrence (using another table similar to ‘probas’). So it only least make the summation.

Step 2: Now that the space is constructed, we need to sum Hamming neighbors, and this is where there’s no shortcut that I know of. You have to assume the space contains randomly generated values. This is the computational bottleneck, and this is where the algorithm I mentioned applies. In fact, as you can see, it doesn’t only go to one space but five, one for each category. So, the sum shouldn’t go to r <= x but should sum the neighbors exactly at x in the corresponding space.

I've edited my code to make the payoff structure dependent on the result of the events, (and therefore it can be dependent on the number of winners for the bet.

It runs slower for now 1900s for 14 matches, of 3 different solutions.

But the structure is just 4 for loops, and the code is quite neat, and I think I kind of see how to apply dynamic programming tricks for memory-speed tradeoff.

Trying to use Pascal's Formula, should probably result in some version of your algorithm.

I'll think about it.

I still have trouble understanding your code. I don’t get why you’re still using ‘probas’ and the inverse probability (line 186:187). I suggest creating an array of 3^14 random numbers and simply using it as input to the function to sum Hamming distances. Since you prefer C++, I made a naive version. I create a space of 3^14 elements, but instead of random numbers, I use sequential ones (the value at position i is equal to i). This makes it easier to compare results. For example, the sum of the Hamming neighbors up to 4 from the first element should give you 18847285404. Let me know how long it takes on your machine (maybe better via GitHub, as it notifies replies to comments). I’m sure you’ll figure out how to optimize it. https://github.com/petopello/PSA/blob/main/naive.cpp

My code is just iterating over the relevant points.

Line 186:187 is the trick to group negative probabilities. You can't apply this if you want to take into account a payout structure which depends on the result.

The more flexible version is Line 421, where you don't use this inverse probability trick, and instead of recomputing proba from scratch could use a precomputed array of 3^14 random element.

Both version should return the same result. You have two boolean useNaive and neighborInBetSpace (4 possibilities) which allows to choose the summation strategy to use, it should give the same results.

I tried understanding your code, but I have trouble understanding the recurrence relation, and where it comes from. In particular how do you make sure you avoid double counting (The neighborhood of a neighborhood contain the original point).

My experiments are about understanding your algorithm in terms of things I understand, namely loop optimizations, partitioning, reordering and memoization. I am trying to make the structure emerge from simple code transformations (to guarantee exactitude). This process usually consist into writing nested loops with everything being recomputed inside the innermost loop, and then find a way to reuse the computations.

Fast code is about memory usage and recomputing from a small cache is often faster than materializing a big array.

To understand recursion, and given that we know the initial value of r (r=4), I suggest replacing the recursive calls with a simple copy-paste of the function body. You’ll see that it’s just four nested loops each one iterating over the 14 dimensions and making changes to them. Each loop starts at the next dimension after the previous one, ensuring that changes are only made in increasing order. This is the key to avoiding duplicates—each pair of modified dimensions is generated only once (e.g., 1&2 is generated, but 2&1 is not). This is controlled by the last parameter (min), which defaults to 0 and sets the starting value of the loops. Another way to see it is that the function returns the sum of neighbors up to distance r, but only by making changes in dimensions equal to or greater than 'min'.

Is the recursion relation you use https://en.wikipedia.org/wiki/Hockey-stick_identity ? Are you sure you can't use https://en.wikipedia.org/wiki/Pascal%27s_rule instead for a simpler recursion. Have you also tried using https://en.wikipedia.org/wiki/Vandermonde%27s_identity ?

Ohh, very interesting, thanks! I’m not entirely sure, but yes, I’d say my naive implementation uses the hockey stick principle. That said, I think many implementations are possible. We could also eliminate recursion by using an index array to store the values of 'min' variable, but I think that would make the behavior harder to understand.

In any case, the key point in my view is that there are 19,321 neighbors at distance ≤4. If we assume as an input condition that their values can be arbitrary—that is, the value of one neighbor has no relation to the others—then regardless of the implementation or mathematical identity used, we’ll end up performing 19,320 summations.

It’s a different story if we want to repeat this process for multiple points. In that case, we can optimize, since some neighbors might be shared and summed only once. This is exactly what my algorithm does: by handling everything in a matrix-based way, it reduces the number of summations per point to just 101 instead of 19,321. I’m not sure if there’s a specific mathematical identity behind this. In fact, I asked on StackExchange but haven’t had much success: https://math.stackexchange.com/questions/5040947/efficient-a...

My algorithm reduces the number of summations by computing all expected values at once, so I can't compute just a single value individually. For a single value, I think the naive approach is the only option. In my first naive implementation, computing 5,000–10,000 values took around 5 minutes, and I didn’t improve it much beyond that.

Once all values are computed, one approach is to search for the highest expected value, but that’s not the only criterion to manage: (1) there's a huge variance issue (you’ll win a lot, but with very low probability), (2) if you're placing multiple bets overlap reduces their combined value.

Great work on your approach! I’ll try to understand the code you linked, but I suspect it’s not doing exactly the same thing (or only for a very specific case). In my code, there's a function called SHN_boxes (which takes ~1s in Colab for this problem), and it’s a shortcut applicable only in some cases (not in this one). Did you use a similar approach?