Training two samples take ~480 seconds

[NAThompson]

Code to reproduce:

#!/usr/bin/env python3
from math import pi as π
import random
import numpy as np
from smt.surrogate_models import RMTB, RMTC
import matplotlib.pyplot as plt


dimension = 4

def generate_sine_surface_data(samples=2):
    training_points = np.full(shape=(samples, dimension), fill_value=np.nan)
    training_values = np.full(shape=samples, fill_value=np.nan)

    for i in range(samples):
        x = np.random.uniform(-1, 1, dimension)
        training_points[i, :] = x
        v = 1
        for j in range(dimension):
            v *= np.sin(π*x[j])
        training_values[i] = v
    return training_points, training_values


def rmtc_surrogate_model(training_points, training_values):
    limits = np.full(shape=(dimension, 2), fill_value=np.nan)
    for i in range(dimension):
        limits[i, 0] = -1
        limits[i, 1] = 1

    model = RMTC(
        xlimits=limits,
        num_elements=10
    )
    model.set_training_values(training_points, training_values)
    model.train()
    return model


if __name__ == "__main__":
    training_points, training_values = generate_sine_surface_data()
    rmtc = rmtc_surrogate_model(training_points, training_values)

Result:

python3 -m cProfile -s time -o rmtc.profile rmtc_issue.py                ✘ 130 571-integer-overflow ◼
___________________________________________________________________________

                                   RMTC
___________________________________________________________________________

 Problem size

      # training points.        : 2

___________________________________________________________________________

 Training

   Training ...
      Pre-computing matrices ...
         Computing dof2coeff ...
         Computing dof2coeff - done. Time (sec): 12.8634088
         Initializing Hessian ...
         Initializing Hessian - done. Time (sec):  0.0033479
         Computing energy terms ...
         Computing energy terms - done. Time (sec): 153.1459041
         Computing approximation terms ...
         Computing approximation terms - done. Time (sec):  3.6191752
      Pre-computing matrices - done. Time (sec): 169.6337998
      Solving for degrees of freedom ...
         Solving initial startup problem (n=234256) ...
            Solving for output 0 ...
               Iteration (num., iy, grad. norm, func.) :   0   0 1.412541449e-01 1.954465085e-02
               Iteration (num., iy, grad. norm, func.) :   0   0 8.538407147e-07 2.031322887e-07
            Solving for output 0 - done. Time (sec): 26.0831141
         Solving initial startup problem (n=234256) - done. Time (sec): 26.0846453
         Solving nonlinear problem (n=234256) ...
            Solving for output 0 ...
               Iteration (num., iy, grad. norm, func.) :   0   0 1.719266581e-06 2.031300540e-07
               Iteration (num., iy, grad. norm, func.) :   0   0 7.548924856e-07 1.614681349e-07
               Iteration (num., iy, grad. norm, func.) :   1   0 4.108689308e-07 1.292814893e-08
               Iteration (num., iy, grad. norm, func.) :   2   0 1.981918673e-07 3.625613917e-09
               Iteration (num., iy, grad. norm, func.) :   3   0 1.380419977e-07 1.567742159e-09
               Iteration (num., iy, grad. norm, func.) :   4   0 3.749522808e-07 1.081064640e-09
               Iteration (num., iy, grad. norm, func.) :   5   0 1.020794481e-07 4.958973451e-10
               Iteration (num., iy, grad. norm, func.) :   6   0 3.519642335e-08 1.290533556e-10
               Iteration (num., iy, grad. norm, func.) :   7   0 2.510705187e-08 2.136727433e-11
               Iteration (num., iy, grad. norm, func.) :   8   0 1.285600568e-08 1.110618394e-11
               Iteration (num., iy, grad. norm, func.) :   9   0 6.092138826e-09 1.248114573e-12
            Solving for output 0 - done. Time (sec): 278.1017668
         Solving nonlinear problem (n=234256) - done. Time (sec): 278.1022830
      Solving for degrees of freedom - done. Time (sec): 304.1938951
   Training - done. Time (sec): 479.8905671

Snakeviz profile screenshotted:

Profile is attached (I had to add a superfluous .txt extension to get github to allow the upload).
rmtc.profile.txt

[relf]

RMTC seems to scale poorly with the num_elements parameter. You could decrease its value and check the quality of the model.

I am not familiar with the RMTC model and its implementation, maybe @hwangjt could give some insights on this topic?

[NAThompson]

Ok, at the suggestion of @relf , I have created a pytest benchmark that reproduces this issue.

First, the here is the runtime as a function of the number of samples; linear scaling:

Now I fix the number of samples at 500, dimension at 3, and then sweep the num_elements:

From this figure it is unclear whether the scaling is quadratic or exponential, but a log-log scale reveals it is a power law:

The scaling with dimension is clearly exponential, and perhaps a little worse than that:

Code to reproduce:

#!/usr/bin/env python3
from math import pi as π
import pathlib
import random
import numpy as np
import pytest
import json
from smt.surrogate_models import RMTB, RMTC
import matplotlib.pyplot as plt


def generate_sine_surface_data(samples=200, dimension=3):
    training_points = np.full(shape=(samples, dimension), fill_value=np.nan)
    training_values = np.full(shape=samples, fill_value=np.nan)

    for i in range(samples):
        x = np.random.uniform(-1, 1, dimension)
        training_points[i, :] = x
        v = 1
        for j in range(dimension):
            v *= np.sin(π*x[j])
        training_values[i] = v
    return training_points, training_values


class Surrogate:
    def __init__(self, samples, num_elements, dimension):
        self.num_elements=num_elements
        self.training_points, self.training_values = generate_sine_surface_data(samples=samples, dimension=dimension)
        self.limits = np.full(shape=(dimension, 2), fill_value=np.nan)
        for i in range(dimension):
            self.limits[i, 0] = -1
            self.limits[i, 1] = 1

    def __call__(self):
        model = RMTC(xlimits=self.limits, num_elements=self.num_elements)
        model.set_training_values(self.training_points, self.training_values)
        model.train()
        return model




@pytest.mark.parametrize('dimension', [1,2,3,4,5])
@pytest.mark.parametrize('num_elements', [4])
@pytest.mark.parametrize('samples', [100])
@pytest.mark.benchmark
def test_num_elements_complexity(benchmark, num_elements, samples, dimension):
    print(f"Dimension = {dimension}")
    s = Surrogate(samples, num_elements, dimension)

    benchmark.pedantic(s, iterations=3, rounds=1)

def runtime_vs_samples(benchmarks):
    samples = []
    runtime = []
    for benchmark in benchmarks:
        params = benchmark['params']
        t = benchmark['stats']['median']
        d = params['dimension']
        num_elements = params['num_elements']
        if num_elements == 3 and d == 3:
            runtime.append(t)
            samples.append(params['samples'])

    
    plt.scatter(samples, runtime)
    plt.xlabel('Samples')
    plt.ylabel('Runtime')
    plt.show()

def runtime_vs_num_elements(benchmarks):
    num_elements = []
    runtime = []
    for benchmark in benchmarks:
        params = benchmark['params']
        t = benchmark['stats']['median']
        d = params['dimension']
        samples = params['samples']
        nelem = params['num_elements']
        runtime.append(t)
        num_elements.append(nelem)

    
    plt.scatter(num_elements, runtime)
    plt.xlabel('num_elements')
    plt.ylabel('Runtime')
    plt.xscale('log')
    plt.yscale('log')
    plt.show()

def runtime_vs_dimension(benchmarks):
    dimension = []
    runtime = []
    for benchmark in benchmarks:
        params = benchmark['params']
        t = benchmark['stats']['median']
        d = params['dimension']
        dimension.append(d)
        runtime.append(t)

    
    plt.scatter(dimension, runtime)
    plt.xlabel('Dimension')
    plt.ylabel('Runtime')
    # plt.xscale('log')
    plt.yscale('log')
    plt.show()



if __name__ == '__main__':
    p = pathlib.Path('~/smt/.benchmarks/Darwin-CPython-3.12-64bit/0017_rmtc.json')
    with open(p) as f:
        j = json.load(f)

    benchmarks = j['benchmarks']
    # runtime_vs_samples(benchmarks)
    # runtime_vs_num_elements(benchmarks)
    runtime_vs_dimension(benchmarks)

Command to reproduce:

$ pip install pytest-benchmark
$ py.test --benchmark-save=rmtc.json rmtc_issue.py

[relf]

Thanks for the thorough analysis!