ENH: LHS based OptimalDesign (scipy.stats.qmc)

[tupui]

This adds a new QMC engine base on optimized LHS.

This was removed from scipy.stats.qmc during its review (#10844).

The optimization loop consists in randomly permuting some points from the design. LHS properties are preserved doing so. All design candidates are compared in terms of centered-C2 discrepancy. This is commonly done in the literature. Also, recomputing the C2 discrepancy is helped with a special method which just update the discrepancy (_perturb_disrepancy from Jin et al. corrected as there where mistakes in the article).

The optimization consists in two loops. An outer loop and an inner loop. For more details, see the reference article: Jin et al., "An efficient algorithm for constructing optimal design of computer experiments." Journal of Statistical Planning and Inference, 2005.

Here is the flowchart of the optimization

[mdhaber]

Does it need to be updated after the LHS bug fix?

[tupui]

Does it need to be updated after the LHS bug fix?

I will double check (and fix the conflict), but I updated it after.

[mdhaber]

Interesting, and seems very useful. I didn't check the implementation of _perturb_discrepancy but it seems to work well in my own test, which I included here. Some documentation suggestions, but probably more important are the questions about the use of basinhopping. Does it perform well on a function that has memory of past inputs? Does "L-BFGS-B" do well with its inputs being rounded?

[mdhaber]

Could use more detail about what func, x0, and bounds are and an example. Example could just be random swaps.

[mdhaber]

Currently, basinhopping does not appear to be performing really well; it seems to explore a very small portion of the possible swaps before terminating. The output design is very similar to the input after one iteration, and they differ only in small-number rows/cols. (e.g. use np.where(before != after)to see the differences).

Whatever optimizer you do end up with, I'd like to see it tested head-to-head against random permutations to see how the discrepancy goes down as a function of time for each. I'm not sure if the overhead of an optimizer is really going to be worth it, since the objective function keeps changing on them. I tested basinhopping against this optimizer:

def swapper(fun, x0, bounds):
    for i in range(1000):
        row_1 = np.random.randint(*bounds[1])
        row_2 = np.random.randint(*bounds[2])
        col = np.random.randint(*bounds[0])
        fun([col, row_1, row_2])
    return x0

and, given the same amount of time, swapper did far better. Here's the (sloppy) script:

import numpy as np
from time import perf_counter_ns
from scipy.stats import qmc

np.random.seed(9392586)
d, n = 5, 400
qmc_gen = qmc.LatinHypercube(d, seed=3295852)
doe = qmc_gen.random(n)
disc = qmc.discrepancy(doe)

allowed_time = 1000000000


print(f"original: {disc}")

t0 = perf_counter_ns()
doe_new = doe.copy()
while(perf_counter_ns() - t0 < allowed_time):
    lhs_gen = qmc.OptimalDesign(d, doe_new, niter=1, seed=20591238)
    doe_new = lhs_gen.random(n)
    disc_new = qmc.discrepancy(doe_new)
    print(f"basinhopping: {disc_new}")

def swapper(fun, x0, bounds):
    for i in range(1000):
        row_1 = np.random.randint(*bounds[1])
        row_2 = np.random.randint(*bounds[2])
        col = np.random.randint(*bounds[0])
        fun([col, row_1, row_2])
    return x0

t0 = perf_counter_ns()
doe_new = doe.copy()
while(perf_counter_ns() - t0 < allowed_time):
    lhs_gen = qmc.OptimalDesign(d, doe_new, method=swapper, niter=1, seed=20591238)
    doe_new = lhs_gen.random(n)
    disc_new = qmc.discrepancy(doe_new)
    print(f"swapper: {disc_new}")

Output:

original: 0.0012317821939307194
basinhopping: 0.001197552675193947
basinhopping: 0.0011950900444781531
basinhopping: 0.0011950900444781531
swapper: 0.0003416202432968696
swapper: 0.00026969876614835187
swapper: 0.00024324343574155805

I wouldn't try super hard to get one of the existing optimizers working. I don't know that any that we have are designed to optimize for a moving target.

Rather than just a finite number of random swaps, a 2-opt like idea is the natural extension, albeit slow.:

def two_opt(fun, x0, bounds):
    thresh = 1e-3  # only start over if improvement is better than thresh
    f0 = fun(x0)
    done = False
    i = 0
    while(not done):
        bound_ranges = [range(bound[0], bound[1]+1) for bound in bounds]
        for col, row_1, row_2 in product(*bound_ranges):
            i += 1
            if not i % 100:  # occasionally print update
                print(col, row_1, row_2)

            f1 = fun([col, row_1, row_2])
            df = f0 - f1
            if df/f0 > thresh:
                print(f"improvement: {df/f0}")
                f0 = f1
                break
        else:
            # thresh /= 10. # could iteratively refine threshold
            # if thresh < 1e-6:
            done = True
        print(f1)
    return x0

[tupui]

This is great. Thanks for working hard on this! I will definitely use your script to test dual_annealing. I will also vary the dimension. Not sure you remember this, but the original code I had for the PR on QMC had 2 strategies. It was first checking the complexity of the problem (just calculating the total number of possible combinations). And if the complexity was bellow a threshold, it was going through all the possible permutations. Otherwise it would use basinhopping.
I am curious to see how dual_annealing performs here. Everyone uses it in the literature so I hope to see better results 🤷‍♂️

[mdhaber]

2-opt is a little different than all possible permutations of the LHS. It is a greedy version of that idea. In the process of making all possible exchanges between two rows for the same column, you don't generate all possible Latin Hypercubes. But yes, best of e.g. 1000 randomly generated Latin Hypercubes and best of all possible Latin Hypercubes (if the problem is small) are also good ideas.

[tupui]

In the first version I also had a mode with a simple loop to generate n LHS and take the best. At the time we discarded this as we said it was too simple and anyone could do this without problem. This could be part of an example to show how discrepancy can be useful.

[tupui]

Interesting metric! I don't have any experience with it, but sounds useful.

I will put this on my list then, thanks.

As for the results, yes, there is some evidence that the match the paper. (Makes me wonder how simple random permutations would fare, but I suppose I can check that on my own when I am looking at this again.)

Did you happen to look into discrepancy as a function of time, or still working on that?

Not yet, just need to find a bit of time for it 😅 But I also would like to be certain that I did not miss anything here.

[tupui]

Ok here is the plot with respect to time, and yes the not fancy method works better (and here I tested lots of combinations of n, dim as I wanted this to work, hum)... I am sending an email to the authors of these papers to see if I have done something wrong. Otherwise I will just go with your idea. Later on we can add a method param if somehow someone finds it 😅

[tupui]

@jschueller @regislebrun if you have some time maybe you can help here 😃 Maybe I have a very stupid bug here like mixing up some variables??

(Just careful about licensing, don't speak about how it's coded in OT)

[tupui]

I just got one reply from G. Damblin and it seems like the algorithm is very difficult to parametrize. He is saying that it should be better in high dimension, but was not that convincing to me. So without further help and dedicated time to research this, I think that this approach will not work for us...

What I propose now is to use your swapper @mdhaber and if someone manages to code ESE and show a real gain (and I am now skeptical about this after having tried... thanks for all your questions which lead me to test this more by the way 😃 ) we can revisit.

(I would not try to push this before the release so we don't need to stress).

[mdhaber]

Sounds good to me. The random row swaps are still quite useful.

[tupui]

Thanks for the review Matt. To summarise:

• I will change the optimizer to dual_annealing. Changed in favour of ESE from the literature.
• Add validation and tests + your extended version for the update.
• Fix small PEP8 and things.
• Do I add a metric for the objective function? Ability to do Linf instead of C2 for instance. Not now.
• Making _perturb_best_doe a parameter? Not so sure as this would require extra logic. Or this could be done using a main optimization class. We would have OptimalDesign(ABC) and OptimalLatinHypercube(OptimalDesign). Does not make sense anymore as very specific algorithm.

The metric is easy and valuable. The other one is more tricky. Optimizing a sample is already something advanced and optimized LHS is commonly used. For the rest it might be too custom so that we have a generic and useful interface.

[mdhaber]

Currently, basinhopping does not appear to be performing really well; it seems to explore a very small portion of the possible swaps before terminating. The output design is very similar to the input after one iteration, and they differ only in small-number rows/cols. (e.g. use np.where(before != after)to see the differences).

Whatever optimizer you do end up with, I'd like to see it tested head-to-head against random permutations to see how the discrepancy goes down as a function of time for each. I'm not sure if the overhead of an optimizer is really going to be worth it, since the objective function keeps changing on them. I tested basinhopping against this optimizer:

def swapper(fun, x0, bounds):
    for i in range(1000):
        row_1 = np.random.randint(*bounds[1])
        row_2 = np.random.randint(*bounds[2])
        col = np.random.randint(*bounds[0])
        fun([col, row_1, row_2])
    return x0

and, given the same amount of time, swapper did far better. Here's the (sloppy) script:

import numpy as np
from time import perf_counter_ns
from scipy.stats import qmc

np.random.seed(9392586)
d, n = 5, 400
qmc_gen = qmc.LatinHypercube(d, seed=3295852)
doe = qmc_gen.random(n)
disc = qmc.discrepancy(doe)

allowed_time = 1000000000


print(f"original: {disc}")

t0 = perf_counter_ns()
doe_new = doe.copy()
while(perf_counter_ns() - t0 < allowed_time):
    lhs_gen = qmc.OptimalDesign(d, doe_new, niter=1, seed=20591238)
    doe_new = lhs_gen.random(n)
    disc_new = qmc.discrepancy(doe_new)
    print(f"basinhopping: {disc_new}")

def swapper(fun, x0, bounds):
    for i in range(1000):
        row_1 = np.random.randint(*bounds[1])
        row_2 = np.random.randint(*bounds[2])
        col = np.random.randint(*bounds[0])
        fun([col, row_1, row_2])
    return x0

t0 = perf_counter_ns()
doe_new = doe.copy()
while(perf_counter_ns() - t0 < allowed_time):
    lhs_gen = qmc.OptimalDesign(d, doe_new, method=swapper, niter=1, seed=20591238)
    doe_new = lhs_gen.random(n)
    disc_new = qmc.discrepancy(doe_new)
    print(f"swapper: {disc_new}")

Output:

original: 0.0012317821939307194
basinhopping: 0.001197552675193947
basinhopping: 0.0011950900444781531
basinhopping: 0.0011950900444781531
swapper: 0.0003416202432968696
swapper: 0.00026969876614835187
swapper: 0.00024324343574155805

I wouldn't try super hard to get one of the existing optimizers working. I don't know that any that we have are designed to optimize for a moving target.

Rather than just a finite number of random swaps, a 2-opt like idea is the natural extension, albeit slow.:

def two_opt(fun, x0, bounds):
    thresh = 1e-3  # only start over if improvement is better than thresh
    f0 = fun(x0)
    done = False
    i = 0
    while(not done):
        bound_ranges = [range(bound[0], bound[1]+1) for bound in bounds]
        for col, row_1, row_2 in product(*bound_ranges):
            i += 1
            if not i % 100:  # occasionally print update
                print(col, row_1, row_2)

            f1 = fun([col, row_1, row_2])
            df = f0 - f1
            if df/f0 > thresh:
                print(f"improvement: {df/f0}")
                f0 = f1
                break
        else:
            # thresh /= 10. # could iteratively refine threshold
            # if thresh < 1e-6:
            done = True
        print(f1)
    return x0

[mdhaber]

I think the interface is a lot better. If users use this feature, they will probably complain and ask for control over the optimization process. I think that's a good thing because it will give us a reason to develop this further. But for now, this is much better than most of the alternatives we've considered because it doesn't change the existing interface in a way that we're unsure about. Adding this one option seems very reasonable to have regardless of where things go, and it's much better than adding an entire class that doesn't really fit the QMCEngine model or even a standalone function which would initially have only one very restricted purpose (to optimize the output of LatinHypercube.random).
Suggested some simplifications to the documentation and code, but looking close.

[mdhaber]

Looks good.

Need to fix the rendering of the new example.

I'd also prefer that the optimization option be case-insensitive. I accidentally tried 'random-cd' and got an error. That annoys me, personally, but I know there are different opinions about that sort of thing.

Here's a simple script to see how the optimization does:

import numpy as np
from scipy.stats import qmc

d, n = 3, 100

rng = np.random.default_rng(9842398122)
lhs_seeds = rng.integers(1000, size=1000)

for seed in lhs_seeds:
    qmc_gen = qmc.LatinHypercube(d, seed=seed)
    doe = qmc_gen.random(n)
    disc = qmc.discrepancy(doe)
    
    qmc_gen = qmc.LatinHypercube(d, seed=seed, optimization='random-CD')
    doe = qmc_gen.random(n)
    optdisc = qmc.discrepancy(doe)
   
    print(optdisc/disc)

It would be nice @tupui if you could produce a few plots showing the discrepancy as a function of iteration so we can gauge whether the convergence criteria are about right. But that is make-better, so we could wait to consider it once we get requests to do so.

But the documentation is OK and the tests seem to check what they need to: the discrepancy update calculation works, the method argument produces a reasonable error message when it's not valid, the optimized design is better than the original, and it's still a Latin Hypercube. So not too much holding this up.

[tupui]

Looks good.

Need to fix the rendering of the new example.

Ooops, should be fixed now.

It would be nice @tupui if you could produce a few plots showing the discrepancy as a function of iteration so we can gauge whether the convergence criteria are about right. But that is make-better, so we could wait to consider it once we get requests to do so.

Ok sure, I will try to make something useful.

But the documentation is OK and the tests seem to check what they need to: the discrepancy update calculation works, the method argument produces a reasonable error message when it's not valid, the optimized design is better than the original, and it's still a Latin Hypercube. So not too much holding this up.

Glad to read, thanks Matt 🥳

[tupui]

@mdhaber I have some results. Seems like I was lucky and 100, 10_000 seems fine.

Note that the values are sometimes not round (ex. 9830) because I am doing the mean over some realizations (here 5 as it does not change much).

d=10, Ns=200

Default values: 10000, 100 | 0.4153428604265345 [  55 9830]
Case 1
1000, 100 | 0.5698480676322689 [   2 1000]
5000, 100 | 0.4499603262234236 [  37 5000]
10000, 100 | 0.4153428604265345 [  55 9830]
20000, 100 | 0.41132065360601955 [  100 10716]
50000, 100 | 0.41132065360601955 [  100 10716]
Case 2
10000, 10 | 0.7494520426748597 [ 10 297]
10000, 100 | 0.4153428604265345 [  55 9830]
10000, 1000 | 0.4144494432252772 [   17 10000]
10000, 2000 | 0.4144494432252772 [   17 10000]
Case 3
50000, 10 | 0.7494520426748597 [ 10 297]
50000, 100 | 0.41132065360601955 [  100 10716]
50000, 1000 | 0.35679136503433057 [   71 50000]
50000, 2000 | 0.35679136503433057 [   71 50000]

-----------------------------------------------------------------

d=5, Ns=100

Default values: 10000, 100 | 0.3119049788602787 [ 100 3947]
Case 1
1000, 100 | 0.3815846896444995 [  10 1000]
5000, 100 | 0.3119049788602787 [ 100 3947]
10000, 100 | 0.3119049788602787 [ 100 3947]
20000, 100 | 0.3119049788602787 [ 100 3947]
50000, 100 | 0.3119049788602787 [ 100 3947]
Case 2
10000, 10 | 0.60966387214425 [ 10 114]
10000, 100 | 0.3119049788602787 [ 100 3947]
10000, 1000 | 0.2778197429019117 [   23 10000]
10000, 2000 | 0.2778197429019117 [   23 10000]
Case 3
50000, 10 | 0.60966387214425 [ 10 114]
50000, 100 | 0.3119049788602787 [ 100 3947]
50000, 1000 | 0.24349852818152812 [  259 48654]
50000, 2000 | 0.24292300734518152 [   72 50000]

-----------------------------------------------------------------

d=2, Ns=100

Default values: 10000, 100 | 0.18105553489676354 [ 100 2330]
Case 1
1000, 100 | 0.19923847532269678 [  33 1000]
5000, 100 | 0.18105553489676354 [ 100 2330]
10000, 100 | 0.18105553489676354 [ 100 2330]
20000, 100 | 0.18105553489676354 [ 100 2330]
50000, 100 | 0.18105553489676354 [ 100 2330]
Case 2
10000, 10 | 0.3887165264896873 [10 74]
10000, 100 | 0.18105553489676354 [ 100 2330]
10000, 1000 | 0.164748928016826 [   97 10000]
10000, 2000 | 0.164748928016826 [   97 10000]
Case 3
50000, 10 | 0.3887165264896873 [10 74]
50000, 100 | 0.18105553489676354 [ 100 2330]
50000, 1000 | 0.16137969129023183 [ 1000 19803]
50000, 2000 | 0.16072167552692784 [ 2000 29274]

Code to reproduce

To reproduce you also need to add self.debug = [n_nochange, n_iters] before _random_cd returns.

import numpy as np
from scipy.stats import qmc

d, n = 10, 200
n_conv = 5

rng = np.random.default_rng(9842398122)
lhs_seeds = rng.integers(1000, size=n_conv)


def opt_disc(seed, n_nochange, n_iters):
    qmc_gen = qmc.LatinHypercube(d, seed=seed)
    doe = qmc_gen.random(n)
    disc = qmc.discrepancy(doe)

    qmc_gen = qmc.LatinHypercube(d, seed=seed, optimization='random-cd')
    qmc_gen._n_nochange = n_nochange
    qmc_gen._n_iters = n_iters
    doe = qmc_gen.random(n)
    optdisc = qmc.discrepancy(doe)

    return np.array([optdisc/disc, *qmc_gen.debug])


n_nochange = 100
n_iters = 10_000

data = [opt_disc(seed, n_nochange=n_nochange, n_iters=n_iters)
         for seed in lhs_seeds]
data = np.mean(data, axis=0)
ratio = data[0]
debug = data[1:].astype(int)
print(f"Default values: {n_iters}, {n_nochange} | {ratio} {debug}")

print("Case 1")

n_nochange = 100
for n_iters in [1000, 5_000, 10_000, 20_000, 50_000]:
    data = [opt_disc(seed, n_nochange=n_nochange, n_iters=n_iters)
            for seed in lhs_seeds]
    data = np.mean(data, axis=0)
    ratio = data[0]
    debug = data[1:].astype(int)
    print(f"{n_iters}, {n_nochange} | {ratio} {debug}")

print("Case 2")
n_iters = 10_000
for n_nochange in [10, 100, 1_000, 2_000]:
    data = [opt_disc(seed, n_nochange=n_nochange, n_iters=n_iters)
            for seed in lhs_seeds]
    data = np.mean(data, axis=0)
    ratio = data[0]
    debug = data[1:].astype(int)
    print(f"{n_iters}, {n_nochange} | {ratio} {debug}")

print("Case 3")
n_iters = 50_000
for n_nochange in [10, 100, 1_000, 2_000]:
    data = [opt_disc(seed, n_nochange=n_nochange, n_iters=n_iters)
            for seed in lhs_seeds]
    data = np.mean(data, axis=0)
    ratio = data[0]
    debug = data[1:].astype(int)
    print(f"{n_iters}, {n_nochange} | {ratio} {debug}")

[mdhaber]

Concerning to use 0d array to avoid an errant type complaint. Is this ever going to return a 0d array?

[tupui]

I tested this and it does return an 0d array if you ask for 0 sample and 0 dim (also in the tests). I can revert if you prefer? (I will overwrite the previous commit if you do)

[mdhaber]

The code

from scipy.stats.qmc import LatinHypercube
lhs = LatinHypercube(0, centered=True)
lhs.random(0)

does not return a 0d array; it returns a 2d array with shape (0, 0). I think that is the correct behavior. I was concerned with the possibility of this function returning a 0d array, which would not be correct.

Still, I am also not sure that we should change the way we write this code to satisfy mypy. It doesn't matter much, but use of a Python float here is what you had in mind originally because it's more natural, and I think the rest of the SciPy codebase agrees.

[tupui]

Sure, agreed. I just reverted this line.

[tupui]

@mdhaber I think I covered all the points. There is just this 0d array question (either the current version is fine or I revert it back to a float) but otherwise all is green now.

[tupui]

In case you did not see @mdhaber, Azure CI is not working right now due to #14529

[mdhaber]

Yes, I noticed. Would you prefer to wait until those checks can run to merge, or do you think that changes since the last successful fun on Azure are unlikely to cause the Azure tests to fail?

[tupui]

Yes, I noticed. Would you prefer to wait until those checks can run to merge, or do you think that changes since the last successful fun on Azure are unlikely to cause the Azure tests to fail?

IMO you can merge. I did not change anything since the last full run on Azure that was failing the type.

scipy/stats/_qmc.py:996: error: Incompatible types in assignment (expression has type "str", variable has type "Optional[Literal['random-cd']]") [assignment]

And the only commit after just add # type: ignore[assignment].

[mdhaber]

Ok! I think this is a valuable change to LatinHypercube. We added the new functionality (optimization of a Latin Hypercube sample to lower the centered discrepancy) with a minimal change to the interface (just the addition of an optional optimization keyword-only argument) to leave open the possibility of expanding in the future. The documentation changes seem complete, and the essential tests have been added:

• accuracy of the method used to update the centered discrepancy w.r.t. elementary permutations of a Latin hypercube
• the optimized Latin hypercube is still a Latin hypercube
• the optimized Latin hypercube has lower centered discrepancy than the original
• basic input validation on the optimization keyword argument

Thanks @tupui!

[tupui]

Thanks Matt for your extensive review which really helped to uncover a lot of things!