Hessian computation/Compatibility with dual numbers

[NAThompson]

Feature Request

Many of the surrogate models provided in the package are $C^2$ and hence their Hessians can in principle be computed. This would be an excellent feature to have; and indeed might be doable via dual numbers.

Of course, dual numbers do not require support on your part; I can simply call your API with them. However, my first attempts to extract Hessians in this way have failed:

#!/usr/bin/env python3
from math import pi as π
import random
import numpy as np
from smt.surrogate_models import RMTB, RMTC
import num_dual # use pip install num_dual


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

    for i in range(samples):
        x = random.uniform(-1, 1)
        y = random.uniform(-2, 2)
        training_points[i, 0] = x
        training_points[i, 1] = y
        training_values[i] = np.sin(π * x) * np.cos(π * y / 2)
    return training_points, training_values


def rmtb_surrogate_model(training_points, training_values):
    limits = np.full(shape=(2, 2), fill_value=np.nan)
    limits[0, 0] = -1
    limits[0, 1] = 1
    limits[1, 0] = -2
    limits[1, 1] = 2

    model = RMTB(
        xlimits=limits,
        order=4,
        num_ctrl_pts=20,
        energy_weight=1e-15,
        regularization_weight=0.0,
        print_global=False,
    )
    model.set_training_values(training_points, training_values)
    model.train()
    return model

if __name__ == "__main__":
    training_points, training_values = generate_sine_surface_data()
    rmtb = rmtb_surrogate_model(training_points, training_values)

    f = lambda x : rmtb.predict_values(np.array([[0.0, 1.0]))[0]
    query_point = np.array([0.3, 0.6])
    out = num_dual.hessian(f, query_point)

Error message:

  File "~/reproduce/lib/python3.11/site-packages/smt/surrogate_models/surrogate_model.py", line 310, in predict_values
    y = self._predict_values(x2)
        ^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/reproduce/lib/python3.11/site-packages/smt/surrogate_models/rmts.py", line 465, in _predict_values
    mtx = self._compute_prediction_mtx(x, 0)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/reproduce/lib/python3.11/site-packages/smt/surrogate_models/rmts.py", line 497, in _compute_prediction_mtx
    data, rows, cols = self._compute_jac_raw(kx, 0, x)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/reproduce/lib/python3.11/site-packages/smt/surrogate_models/rmtb.py", line 104, in _compute_jac_raw
    t[:, kx] = (x[:, kx] - xlimits[kx, 0]) / (xlimits[kx, 1] - xlimits[kx, 0])
    ~^^^^^^^
TypeError: float() argument must be a string or a real number, not 'builtins.Dual2Vec64'

Can Hessians be supported directly?

[relf]

Thanks, I was not aware of the num_dual library. I can not run the code above with the expected error. What is the type/value of query_points[0]?

Anyway seems it is not compatible with num_dual! Regarding hessians implementations, I would say we are always open to contributions.

[NAThompson]

I can not run the code above with the expected error. What is the type/value of query_points[0]?

Sorry-I introduced a typo in the reproducer while I was trying to pare it down. Fixed now.

[relf]

I had to change the code of f to (try to) make it work (missing ] , try to use item() to get a "scalar"):

    def f(x):
        y = rmtb.predict_values(np.array([[0.0, 1.0]])).item()
        print(y, type(y))
        return y

When I run it:

> python test_num_dual.py
-3.949565971416e-06 <class 'float'>
Traceback (most recent call last):
  File "D:\rlafage\test_num_dual.py", line 52, in <module>
    out = num_dual.hessian(f, query_point)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: argument 'f' must return a scalar.

I do not understand this error and I do not get yours, anyway still not ready to use num_dual 😉

[NAThompson]

That looks like it may well be a bug in num_dual . . . you've made a pretty compelling case that f returns a scalar . . .