[Feature Request] Hypervolume computation with BoTorch is thousands of times slower than with moocore

[MLopez-Ibanez]

🚀 Feature Request

Motivation

https://botorch.org/api/utils.html#botorch.utils.multi_objective.hypervolume.Hypervolume

says: TODO: write this in C++ for faster looping.

There is a C (C++ compatible) implementation available here: https://github.com/MLopez-Ibanez/eaf/blob/master/src/mo-tools/hv.c

[sdaulton]

Yeah, that would be nice. Using box decompositions (DominatedPartitioning) to compute hypervolumes seems to work more efficiently than the dimension-sweep algorithm, empirically. Both would likely benefit from C++

[Balandat]

We'd also need this to be auto-differentiable in pytorch though. So the idea was to (maybe) write this using the pytorch C++ frontend as an extension.

[sdaulton]

The current implementation is not differentiable and is simply used for evaluation purposes.

[Balandat]

Oh I see different implementation, my bad.

[MLopez-Ibanez]

Yeah, that would be nice. Using box decompositions (DominatedPartitioning) to compute hypervolumes seems to work more efficiently than the dimension-sweep algorithm, empirically. Both would likely benefit from C++

Maybe. I believe the state-of-the-art is pretty well summarised here: https://dl.acm.org/doi/fullHtml/10.1145/3453474#tab1 And there it recommends the 2D and 3D dimension-sweep algorithms for 2D and 3D. The box-decomposition is only recommended for 5D and 6D.

In any case, I would be surprised if the Python implementation of either algorithm can get anywhere close to the C version of dimension-sweep for 2D and 3D.

It would be nice to have a highly-optimized C/C++ library of fundamental algorithms (not only hypervolume, but dominance filtering, nondominated sorting, other metrics like IGD+ and epsilon, etc) for multi-objective optimization that can be used in R and Python.

[MLopez-Ibanez]

There is a C (C++ compatible) implementation available here: https://github.com/MLopez-Ibanez/eaf/blob/master/src/mo-tools/hv.c

By the way, I'm happy to modify the C code above to make it easier to integrate into BoTorch.

[Balandat]

It would be nice to have a highly-optimized C/C++ library of fundamental algorithms (not only hypervolume, but dominance filtering, nondominated sorting, other metrics like IGD+ and epsilon, etc) for multi-objective optimization that can be used in R and Python.

Indeed, that would be nice.

By the way, I'm happy to modify the C code above to make it easier to integrate into BoTorch.

BoTorch currently doesn't contain any compiled code, which simplifies building and packaging a lot, and so I'd only want to add this to BoTorch if it brings substantial benefits. That said, if this lived in a separate library that is well implemented and maintained I'd be happy to add it as a (potentially optional) dependency. This might more sense especially if the scope for this work is broader and potentially includes more general python and R bindings.

[MLopez-Ibanez]

I did some quick experiments comparing the BoTorch implementation (botorch-0.12.0) of the hypervolume to the one provided by moocore (https://multi-objective.github.io/moocore/python/), which is implemented in C.

It seems BoTorch is several thousands times slower. Maybe I am doing something wrong? I used this code.

import numpy as np
import moocore

from botorch.utils.multi_objective.hypervolume import Hypervolume as botorch_HV
import torch

import timeit

def read_data(name):
    files = {"DTLZLinearShape.3d": "~/work/perfassess/hyperv/uc/test/inp/DTLZLinearShape.3d.front.1000pts.10",
             "DTLZLinearShape.4d" : "~/work/perfassess/hyperv/uc/test/inp/DTLZLinearShape.4d.front.1000pts.10"}
    x = moocore.read_datasets(files[name])[:, :-1]
    x = moocore.filter_dominated(x, maximise=True)
    return x


name = "DTLZLinearShape.3d"
x = read_data(name)

ref = np.zeros(x.shape[1])
botorch_hv = botorch_HV(ref_point=torch.from_numpy(ref))
moocore_hv = moocore.Hypervolume(ref = ref, maximise = True)

n = np.arange(500, 3000 + 1, 500)

botorch_values = []
moocore_values =  []
botorch_times = []
moocore_times =  []
for maxrow in n:
    z = x[:maxrow, :]
    z_torch = torch.from_numpy(z)
    botorch_values += [botorch_hv.compute(z_torch)]
    moocore_values += [moocore_hv(z)]
    duration = timeit.timeit('botorch_hv.compute(z_torch)', globals=globals(), number = 2)  
    botorch_times += [duration]
    duration = timeit.timeit('moocore_hv(z)', globals=globals(), number = 2)
    moocore_times += [duration]
    
botorch_values = np.array(botorch_values)
moocore_values = np.array(moocore_values)
assert np.allclose(botorch_values, moocore_values)

botorch_times = np.array(botorch_times)
moocore_times = np.array(moocore_times)

import pandas as pd
import matplotlib.pyplot as plt
cpu = "Intel i5-6200U 2.30GHz"
df = pd.DataFrame(dict(n = n, Ratio = botorch_times/moocore_times)).set_index('n')
df.plot(grid=True, ylabel='Botorch/moocore (seconds)', title = f'HV computation for {name} ({cpu})')
plt.show()

[eytan]

I'll let somebody else weigh in on whether this is the right way to benchmark, but for reference, which algorithm are you using for computing HV in moocore?

…

On Thu, Nov 28, 2024 at 12:34 PM Manuel López-Ibáñez < ***@***.***> wrote: I did some quick experiments comparing the BoTorch implementation (botorch-0.12.0) of the hypervolume to the one provided by moocore ( https://multi-objective.github.io/moocore/python/), which is implemented in C. Figure_1.png (view on web) <https://github.com/user-attachments/assets/be4d224e-025c-4382-85a0-c97928615b94> It seems BoTorch is several thousands times slower. Maybe I am doing something wrong? I used this code. import numpy as npimport moocore from botorch.utils.multi_objective.hypervolume import Hypervolume as botorch_HVimport torch import timeit def read_data(name): files = {"DTLZLinearShape.3d": "~/work/perfassess/hyperv/uc/test/inp/DTLZLinearShape.3d.front.1000pts.10", "DTLZLinearShape.4d" : "~/work/perfassess/hyperv/uc/test/inp/DTLZLinearShape.4d.front.1000pts.10"} x = moocore.read_datasets(files[name])[:, :-1] x = moocore.filter_dominated(x, maximise=True) return x name = "DTLZLinearShape.3d"x = read_data(name) ref = np.zeros(x.shape[1])botorch_hv = botorch_HV(ref_point=torch.from_numpy(ref))moocore_hv = moocore.Hypervolume(ref = ref, maximise = True) n = np.arange(500, 3000 + 1, 500) botorch_values = []moocore_values = []botorch_times = []moocore_times = []for maxrow in n: z = x[:maxrow, :] z_torch = torch.from_numpy(z) botorch_values += [botorch_hv.compute(z_torch)] moocore_values += [moocore_hv(z)] duration = timeit.timeit('botorch_hv.compute(z_torch)', globals=globals(), number = 2) botorch_times += [duration] duration = timeit.timeit('moocore_hv(z)', globals=globals(), number = 2) moocore_times += [duration] botorch_values = np.array(botorch_values)moocore_values = np.array(moocore_values)assert np.allclose(botorch_values, moocore_values) botorch_times = np.array(botorch_times)moocore_times = np.array(moocore_times) import pandas as pdimport matplotlib.pyplot as pltcpu = "Intel i5-6200U 2.30GHz"df = pd.DataFrame(dict(n = n, Ratio = botorch_times/moocore_times)).set_index('n')df.plot(grid=True, ylabel='Botorch/moocore (seconds)', title = f'HV computation for {name} ({cpu})')plt.show() — Reply to this email directly, view it on GitHub <#1685 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAAW34KYSH52YGJX6OHB6D32C5HZFAVCNFSM6AAAAABSVTGDASVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKMBWGU3DENJUHE> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[MLopez-Ibanez]

Documented here: https://multi-objective.github.io/moocore/python/reference/generated/moocore.hypervolume.html#moocore.hypervolume

[Balandat]

I haven't looked into the details, but I wouldn't be surprised if the implementation of this particular algorithm were indeed many orders of magnitude slower. This algorithm is very loopy so this is going to be super slow in python to begin with. On top of that this is done in pytorch so there is probably a lot of overhead from building the (very complex) compute graph in the background. This algorithm in botorch really shouldn't be used (and it's not used in our multi-objective acquisition functions). We should probably make this more clear in the docs and potentially just remove this algorithm.

[MLopez-Ibanez]

It would be nice to have a highly-optimized C/C++ library of fundamental algorithms (not only hypervolume, but dominance filtering, nondominated sorting, other metrics like IGD+ and epsilon, etc) for multi-objective optimization that can be used in R and Python.

Indeed, that would be nice.

Well, that is what moocore is meant to be.

[MLopez-Ibanez]

https://botorch.org/api/utils.html#botorch.utils.multi_objective.pareto.is_non_dominated

versus

https://multi-objective.github.io/moocore/python/reference/generated/moocore.is_nondominated.html

import numpy as np
import moocore
import torch
import timeit

from botorch.utils.multi_objective.pareto import is_non_dominated as botorch_is_non_dominated

name = "CPFs.txt"
x = moocore.get_dataset(name)[:, :-1]
botorch_times = []
moocore_times =  []
n = np.arange(100, len(x) + 1, 250)
botorch_times = []
moocore_times =  []
for maxrow in n:
    z = x[:maxrow, :]
    z_torch = torch.from_numpy(z)
    botorch_values = botorch_is_non_dominated(z_torch)
    duration = timeit.timeit('botorch_is_non_dominated(z_torch)', globals=globals(), number = 5)
    botorch_times += [duration]
    moocore_values = moocore.is_nondominated(z, maximise=True)
    duration = timeit.timeit('moocore.is_nondominated(z, maximise=True)', globals=globals(), number = 5)
    moocore_times += [duration]
    assert np.allclose(botorch_values,moocore_values)

botorch_times = np.array(botorch_times)
moocore_times = np.array(moocore_times)

import pandas as pd
import matplotlib.pyplot as plt
cpu = "Intel i5-6200U 2.30GHz"
df = pd.DataFrame(dict(n = n, BoTorch = botorch_times, moocore = moocore_times)).set_index('n')
df.plot(grid=True, ylabel='CPU time (seconds)', title = f'is_non_dominated() for {name} ({cpu})')
plt.yscale("log")
plt.show()