Module func: Functional TT-format including Chebyshev interpolation

Package teneva, module func: functional TT-format with interpolation.

This module contains the functions for construction of the functional TT-representation, including Chebyshev interpolation in the TT-format, as well as calculating the values of the function using the constructed interpolation coefficients.



teneva.func.func_basis(X, m=10, kind='cheb')[source]

Compute the basis functions in the given points.

The function computes values of the first m basis functions (Chebyshev polynomials in the current version) in the given points X.

Parameters:

  • X (np.ndarray) – spatial points of interest (it is 2D array of the shape [samples, d], where d is a number of dimensions). All points should be from the interval [-1, 1] for “cheb” kind and from the interval [0, pi] for “sin” kind.

  • m (int) – maximum order of the basis function (>= 1). The first m functions (of the order 0, 1, …, m-1) will be computed.

  • kind (str) – kind of the basis (“cheb”).

Returns:

values of the basis functions of the order 0, 1, …, m-1 in X points (it is 3D array of the shape [m, samples, d]).

Return type:

np.ndarray

Examples:

X = np.array([           # Two 4-dim points
    [0., 0., 0., 0.],
    [1., 1., 1., 1.],
])

m = 3                    # Maximum order of polynomial

# Compute Chebyshev polynomials:
T = teneva.func_basis(X, m)

print(T.shape)
print(T)

# >>> ----------------------------------------
# >>> Output:

# (3, 2, 4)
# [[[ 1.  1.  1.  1.]
#   [ 1.  1.  1.  1.]]
#
#  [[ 0.  0.  0.  0.]
#   [ 1.  1.  1.  1.]]
#
#  [[-1. -1. -1. -1.]
#   [ 1.  1.  1.  1.]]]
#


teneva.func.func_diff_matrix(a, b, n, m=1, kind='cheb')[source]

Construct the differential matrix (Chebyshev or Sin) of any order.

The function returns the matrix D (if m=1), which, for the known vector y of values of a one-dimensional function on the related grid (“Chebyshev” grid if kind is “cheb” or uniform grid if kind is “sin”), gives its first derivative, i.e., y’ = D y. If the argument m is greater than 1, then the function returns a list of matrices corresponding to the first derivative, the second derivative, and so on. Note that the derivative error can be significant near the boundary of the region.

Parameters:

  • a (float) – grid lower bound.

  • b (float) – grid upper bound.

  • n (int) – grid size.

  • m (int) – the maximum order of derivative.

  • kind (str) – kind of the basis (“cheb” or “sin”).

Returns:

the differential matrices of order 1, 2, …, m if m > 1, or only one matrix corresponding to the first derivative if m = 1.

Return type:

list of np.ndarray or np.ndarray

Examples:

Let build an analytic function for demonstration:

a = -2.   # Grid lower bound
b = +3.   # Grid upper bound
n = 1000  # Grid size

# Function and its first derivative:
f     = lambda x: np.sin(x**3) + np.exp(-x**2)
f_der = lambda x: 3. * x**2 * np.cos(x**3) - 2. * x * np.exp(-x**2)

# Chebyshev grid and function values on the grid:
i = np.arange(n)
x = teneva.ind_to_poi(i, a, b, n, kind='cheb')
y = f(x)

We can compute the derivative for “y” by Chebyshev differential matrix:

D = teneva.func_diff_matrix(a, b, n)
z = D @ y

Let check the result:

z_real = f_der(x)

e_nrm = np.linalg.norm(z - z_real) / np.linalg.norm(z_real)
e_max = np.max(np.abs((z - z_real) / z_real))

print(f'Error nrm : {e_nrm:-7.1e}')
print(f'Error max : {e_max:-7.1e}')

# >>> ----------------------------------------
# >>> Output:

# Error nrm : 7.5e-13
# Error max : 1.4e-09
#

We can also calculate higher order derivatives:

D1, D2, D3 = teneva.func_diff_matrix(a, b, n, m=3)
z = [D1 @ y, D2 @ y, D3 @ y]

Let check the result:

z1_real = 3. * x**2 * np.cos(x**3) - 2. * x * np.exp(-x**2)

z2_real = 6. * x * np.cos(x**3) - 9. * x**4 * np.sin(x**3)
z2_real += - 2. * np.exp(-x**2) + 4. * x**2 * np.exp(-x**2)

z3_real = 6. * np.cos(x**3) - 18. * x**3 * np.sin(x**3)
z3_real += - 36. * x**3 * np.sin(x**3) - 27. * x**6 * np.cos(x**3)
z3_real += 4. * x * np.exp(-x**2)
z3_real += 8. * x * np.exp(-x**2) - 8. * x**3 * np.exp(-x**2)

z_real = [z1_real, z2_real, z3_real]

for k in range(3):
    e_nrm = np.linalg.norm(z[k] - z_real[k]) / np.linalg.norm(z_real[k])
    e_max = np.max(np.abs((z[k] - z_real[k]) / z_real[k]))
    print(f'Der # {k+1} | Error nrm : {e_nrm:-7.1e} | Error max : {e_max:-7.1e}')

# >>> ----------------------------------------
# >>> Output:

# Der # 1 | Error nrm : 7.5e-13 | Error max : 1.4e-09
# Der # 2 | Error nrm : 4.9e-09 | Error max : 4.3e-08
# Der # 3 | Error nrm : 1.3e-05 | Error max : 1.4e-03
#

We may also use the “sin” basis (DRAFT!!!):

a = 0.       # Grid lower bound
b = np.pi    # Grid upper bound
n = 1000     # Grid size

# Function and its first derivative:
f     = lambda x: np.sin(x)
f_der = lambda x: np.cos(x)

# Uniform grid and function values on the grid:
i = np.arange(n)
x = teneva.ind_to_poi(i, a, b, n, kind='uni')
y = f(x)

D = teneva.func_diff_matrix(a, b, n, kind='sin')
z = D @ y

z_real = f_der(x)

e_nrm = np.linalg.norm(z - z_real) / np.linalg.norm(z_real)
e_max = np.max(np.abs((z - z_real) / z_real))

print(f'Error nrm : {e_nrm:-7.1e}')
print(f'Error max : {e_max:-7.1e}')

# >>> ----------------------------------------
# >>> Output:

# Error nrm : 5.3e+02
# Error max : 3.2e+05
#


teneva.func.func_get(X, A, a=None, b=None, z=0.0, funcs=None, kind='cheb', skip_out=None)[source]

Compute the functional TT-approximation in given points (approx. f(X)).

Parameters:

  • X (np.ndarray) – spatial points of interest (it is 2D array of the shape [samples, d], where d is the number of dimensions). In the case of only one point, it may be also 1D array.

  • A (list) – TT-tensor of the interpolation coefficients (it has d dimensions).

  • a (float, list, np.ndarray) – grid lower bounds for each dimension (list or np.ndarray of length d). It may be also float, then the lower bounds for each dimension will be the same.

  • b (float, list, np.ndarray) – grid upper bounds for each dimension (list or np.ndarray of length d). It may be also float, then the upper bounds for each dimension will be the same.

  • z (float) – the value for points, which are outside the spatial grid.

  • funcs (list of callables) – optional custom basis functions.

  • kind (str) – kind of the basis (“cheb”). It is not used, kept for down compatibility.

  • skip_out (bool) – if flag is set, then the values outside the spatial grid will be set to z values.

Returns:

approximated function values in X points (it is 1D array of the shape [samples]). In the case of only one input point X, it will be float value.

Return type:

np.ndarray

Examples:

# In the beginning we compute the function values on the Chebyshev
# grid using TT-cross (see cross function for more details):
from scipy.optimize import rosen
f = lambda X: rosen(X.T)        # Target function

a = [-2., -4., -3., -2.]        # Grid lower bounds
b = [+2., +3., +4., +2.]        # Grid upper bounds
n = [5, 6, 7, 8]                # Grid size
Y0 = teneva.rand(n, r=2)        # Initial approximation for TT-cross
e = 1.E-3                       # Accuracy for TT-CROSS
eps = 1.E-6                     # Accuracy for truncation

Y = teneva.cross(lambda I: f(teneva.ind_to_poi(I, a, b, n, 'cheb')),
    Y0, e=e, log=True)
Y = teneva.truncate(Y, eps)
teneva.show(Y)

# >>> ----------------------------------------
# >>> Output:

# # pre | time:      0.003 | evals: 0.00e+00 | rank:   2.0 |
# #   1 | time:      0.010 | evals: 3.12e+02 | rank:   4.0 | e: 7.9e+04 |
# #   2 | time:      0.015 | evals: 1.09e+03 | rank:   6.0 | e: 0.0e+00 | stop: e |
# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

# Then we should compute the TT-tensor for Chebyshev interpolation
# coefficients (see func_int function for more details):
A = teneva.func_int(Y)

teneva.show(A) # Show the result

# >>> ----------------------------------------
# >>> Output:

# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

# Finally we compute the approximation in selected points inside
# the bounds (the values for points outside the bounds will be set as "z"):
X = np.array([
    [0., 0., 0., 0.],
    [0., 2., 3., 2.],
    [1., 1., 1., 1.],
    [1., 1., 1., 99999999],
])

Z = teneva.func_get(X, A, a, b, z=-1.)

print(Z)    # Print the result
print(f(X)) # We can check the result by comparing it to the true values

# >>> ----------------------------------------
# >>> Output:

# [ 3.00000000e+00  5.40600000e+03  7.73070497e-12 -1.00000000e+00]
# [3.0000000e+00 5.4060000e+03 0.0000000e+00 9.9999996e+17]
#

We may also compute the value for only one point:

x = np.array([0., 2., 3., 2.])
z = teneva.func_get(x, A, a, b, z=-1.)

print(z)    # Print the result
print(f(x)) # We can check the result by comparing it to the true value

# >>> ----------------------------------------
# >>> Output:

# 5405.999999999997
# 5406.0
#


teneva.func.func_gets(A, m=None, kind='cheb')[source]

Compute the functional approximation (TT-tensor) all over the new grid.

Parameters:

  • A (list) – TT-tensor of the interpolation coefficients (it has d dimensions).

  • m (int, float, list, np.ndarray) – tensor size for each dimension of the new grid (list or np.ndarray of length d). It may be also int/float, then the size for each dimension will be the same. If it is not set, then original grid size (from the interpolation) will be used.

  • kind (str) – kind of the basis (“cheb” or “sin”).

Returns:

TT-tensor of the approximated function values on the full new grid. This relates to the d-dimensional array of the shape m.

Return type:

list

Note

Sometimes additional rounding of the result is relevant. Use for this Z = teneva.truncate(Z, e) (e.g., e = 1.E-8) after the function call.

Examples:

# In the beginning we compute the function values on the Chebyshev
# grid using TT-cross (see cross function for more details):
from scipy.optimize import rosen
f = lambda X: rosen(X.T)        # Target function

a = [-2., -4., -3., -2.]        # Grid lower bounds
b = [+2., +3., +4., +2.]        # Grid upper bounds
n = [5, 6, 7, 8]                # Grid size
Y0 = teneva.rand(n, r=2)        # Initial approximation for TT-cross
e = 1.E-3                       # Accuracy for TT-CROSS
eps = 1.E-6                     # Accuracy for truncation


Y = teneva.cross(lambda I: f(teneva.ind_to_poi(I, a, b, n, 'cheb')),
    Y0, e=e, log=True)
Y = teneva.truncate(Y, eps)
teneva.show(Y)

# >>> ----------------------------------------
# >>> Output:

# # pre | time:      0.005 | evals: 0.00e+00 | rank:   2.0 |
# #   1 | time:      0.016 | evals: 3.12e+02 | rank:   4.0 | e: 6.7e+04 |
# #   2 | time:      0.028 | evals: 1.09e+03 | rank:   6.0 | e: 4.3e-09 | stop: e |
# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

# Then we should compute the TT-tensor for Chebyshev interpolation
# coefficients (see func_int function for more details):
A = teneva.func_int(Y)

teneva.show(A) # Show the result

# >>> ----------------------------------------
# >>> Output:

# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

m = [7, 8, 9, 10] # New size of the grid

# Compute tensor on finer grid:
Z = teneva.func_gets(A, m)

teneva.show(Z)

# >>> ----------------------------------------
# >>> Output:

# TT-tensor     4D : |7| |8| |9| |10|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

# We can compute interpolation coefficients on the new grid:
B = teneva.func_int(Z)

teneva.show(B)

# >>> ----------------------------------------
# >>> Output:

# TT-tensor     4D : |7| |8| |9| |10|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

# Finally we compute the approximation in selected points
# inside the bounds for 2 different approximations:
X = np.array([
    [0., 0., 0., 0.],
    [0., 2., 3., 2.],
    [1., 1., 1., 1.],
    [1., 1., 1., 99999999],
])

z1 = teneva.func_get(X, A, a, b, z=-1.)
z2 = teneva.func_get(X, B, a, b, z=-1.)

# We can check the result by comparing it to the true values:
print(z1)
print(z2)
print(f(X))

# >>> ----------------------------------------
# >>> Output:

# [ 3.00000000e+00  5.40600000e+03  7.73070497e-12 -1.00000000e+00]
# [ 3.00000000e+00  5.40600000e+03  1.20365939e-11 -1.00000000e+00]
# [3.0000000e+00 5.4060000e+03 0.0000000e+00 9.9999996e+17]
#

We may also use “sin” basis:

Y = teneva.cross(lambda I: f(teneva.ind_to_poi(I, a, b, n, 'uni')),
    Y0, e=e, log=True)
Y = teneva.truncate(Y, eps)
A = teneva.func_int(Y, kind='sin')
Z = teneva.func_gets(A, m, kind='sin')
teneva.show(Z)

# >>> ----------------------------------------
# >>> Output:

# # pre | time:      0.004 | evals: 0.00e+00 | rank:   2.0 |
# #   1 | time:      0.013 | evals: 3.12e+02 | rank:   4.0 | e: 5.7e+04 |
# #   2 | time:      0.019 | evals: 1.09e+03 | rank:   6.0 | e: 0.0e+00 | stop: e |
# TT-tensor     4D : |7| |8| |9| |10|
# <rank>  =    3.0 :   \3/ \3/ \3/
#


teneva.func.func_int(Y, kind='cheb')[source]

Compute the TT-tensor for functional interpolation coefficients.

Parameters:

Y (list) – TT-tensor with function values on the corresponding grid.

Returns:

TT-tensor that collects interpolation coefficients. It has the same shape as the given tensor Y.

Return type:

list

Note

Sometimes additional rounding of the result is relevant. Use for this A = teneva.truncate(A, e) (e.g., e = 1.E-8) after the function call.

Examples:

# In the beginning we compute the function values on the Chebyshev
# grid using TT-cross (see cross function for more details):
from scipy.optimize import rosen
f = lambda X: rosen(X.T)        # Target function

a = [-2., -4., -3., -2.]        # Grid lower bounds
b = [+2., +3., +4., +2.]        # Grid upper bounds
n = [5, 6, 7, 8]                # Grid size
Y0 = teneva.rand(n, r=2)        # Initial approximation for TT-cross
e = 1.E-3                       # Accuracy for TT-CROSS
eps = 1.E-6                     # Accuracy for truncation


Y = teneva.cross(lambda I: f(teneva.ind_to_poi(I, a, b, n, 'cheb')),
    Y0, e=e, log=True)
Y = teneva.truncate(Y, eps)
teneva.show(Y)

# >>> ----------------------------------------
# >>> Output:

# # pre | time:      0.004 | evals: 0.00e+00 | rank:   2.0 |
# #   1 | time:      0.010 | evals: 3.12e+02 | rank:   4.0 | e: 1.1e+05 |
# #   2 | time:      0.016 | evals: 1.09e+03 | rank:   6.0 | e: 6.0e-09 | stop: e |
# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

# Then we can compute the TT-tensor for Chebyshev
# interpolation coefficients:
A = teneva.func_int(Y)

teneva.show(A)

# >>> ----------------------------------------
# >>> Output:

# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

We may also use “sin” basis:

Y = teneva.cross(lambda I: f(teneva.ind_to_poi(I, a, b, n, 'uni')),
    Y0, e=e, log=True)
Y = teneva.truncate(Y, eps)
teneva.show(Y)

# >>> ----------------------------------------
# >>> Output:

# # pre | time:      0.004 | evals: 0.00e+00 | rank:   2.0 |
# #   1 | time:      0.010 | evals: 3.12e+02 | rank:   4.0 | e: 9.5e+04 |
# #   2 | time:      0.018 | evals: 1.09e+03 | rank:   6.0 | e: 1.4e-08 | stop: e |
# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#

# Then we can compute the TT-tensor for Sin
# interpolation coefficients:
A = teneva.func_int(Y, kind='sin')

teneva.show(A)

# >>> ----------------------------------------
# >>> Output:

# TT-tensor     4D : |5| |6| |7| |8|
# <rank>  =    3.0 :   \3/ \3/ \3/
#


teneva.func.func_int_general(Y, X, basis_func, rcond=1e-06)[source]

Compute the TT-tensor for functional interpolation coefficients.

Parameters:

  • Y (list) – TT-tensor.

  • X (list, np.ndarray) – values of continuous argument for each TT-core. It should be 2-dim array or 1-dim array if the values are the same for all TT-cores.

  • basis_func (function) – function, which corresponds to the values of the basis functions in the TT-format. It should return np.ndarray of the size n x m, where n is a number of one-dimensional points and m is a number of basis functions.

Returns:

TT-tensor, which represents the interpolation coefficients in terms of the functional TT-format.

Return type:

list

Examples:

a = -2. # Lower bound for continuous grid
b = +3. # Upper bound for continuous grid
d = 4   # Dimension of the grid
n = 10  # Number of grid points

# Build grid points:
I = np.arange(n)
X = teneva.ind_to_poi(I, a, b, n)

# Random TT-tensor:
Y = teneva.rand([n]*d, r=4)

# basis_func = TODO
# A = teneva.func_int_general(Y, X, basis_func, rcond=1e-6)


teneva.func.func_sum(A, a, b, kind='cheb')[source]

Integrate a function from its functional approximation in the TT-format.

Parameters:

  • A (list) – TT-tensor of the interpolation coefficients (it has d dimensions).

  • a (float, list, np.ndarray) – grid lower bounds for each dimension (list or np.ndarray of length d). It may be also float, then the lower bounds for each dimension will be the same.

  • b (float, list, np.ndarray) – grid upper bounds for each dimension (list or np.ndarray of length d). It may be also float, then the upper bounds for each dimension will be the same.

  • kind (str) – kind of the basis (“cheb” or “sin”).

Returns:

the value of the integral.

Return type:

float

Note

This function works only for symmetric grids!

Examples:

# In the beginning we compute the function values on the Chebyshev
# grid using TT-cross (see cheb_bld function for more details):

d = 4
def f(X): # Target function
    a = 2.
    r = np.exp(-np.sum(X*X, axis=1) / a) / (np.pi * a)**(d/2)
    return r.reshape(-1)

a = [-12., -14., -13., -11.]    # Grid lower bounds
b = [+12., +14., +13., +11.]    # Grid upper bounds
n = [50, 50, 50, 50]            # Grid size
Y0 = teneva.rand(n, r=2)        # Initial approximation for TT-cross
e = 1.E-5                       # Accuracy for TT-CROSS
eps = 1.E-6                     # Accuracy for truncation

Y = teneva.cross(lambda I: f(teneva.ind_to_poi(I, a, b, n, 'cheb')),
    Y0, e=e, log=True)
Y = teneva.truncate(Y, eps)
teneva.show(Y)

# >>> ----------------------------------------
# >>> Output:

# # pre | time:      0.004 | evals: 0.00e+00 | rank:   2.0 |
# #   1 | time:      0.013 | evals: 2.40e+03 | rank:   4.0 | e: 1.0e+00 |
# #   2 | time:      0.029 | evals: 8.40e+03 | rank:   6.0 | e: 7.4e-09 | stop: e |
# TT-tensor     4D : |50| |50| |50| |50|
# <rank>  =    1.0 :    \1/  \1/  \1/
#

# Then we should compute the TT-tensor for Chebyshev interpolation
# coefficients (see func_int function for more details):
A = teneva.func_int(Y)

teneva.show(A)

# >>> ----------------------------------------
# >>> Output:

# TT-tensor     4D : |50| |50| |50| |50|
# <rank>  =    1.0 :    \1/  \1/  \1/
#

# Finally we compute the integral:
v = teneva.func_sum(A, a, b)

print(v) # Print the result (the real value is 1.)

# >>> ----------------------------------------
# >>> Output:

# 1.000000019159871
#