imate.InterpolateLogdet#

class imate.InterpolateLogdet(A, B=None, options={}, verbose=False, kind='imbf', ti=[], **kwargs)#

Interpolate the log-deterinant of an affine matrix function.

Parameters:

Anumpy.ndarray, scipy.sparse matrix

Symmetric positive-definite matrix (positive-definite if p is non-positive). Matrix can be dense or sparse.

Warning

Symmetry and positive (semi-) definiteness of A will not be checked. Make sure A satisfies these conditions.

Bnumpy.ndarray, scipy.sparse matrix, default=None

Symmetric positive-definite matrix (positive-definite if p is non-positive). Matrix can be dense or sparse. B should have the same size and type of A. If B is None (default value), it is assumed that B is the identity matrix.

Warning

Symmetry and positive (semi-) definiteness of B will not be checked. Make sure B satisfies these conditions.

optionsdict, default={}

At each interpolation point \(t_i\), the logdeterminant is computed using the function imate.logdet(). The options passes a dictionary of arguments to this functions.

verbosebool, default=False

If True, it prints some information about the computation process.

kind{‘ext’, ‘eig’, ‘mbf’, ‘imbf’, ‘rbf’, ‘crf’, ‘spl’, ‘rpf’}, default: ‘imbf’

The algorithm of interpolation, which are the same as the algorithms in imate.InterpolateSchatten. See documentation for each specific algorithm:

tifloat or array_like(float), default=None

Interpolation points, which can be a single point, a list, or an array of points. The interpolator honors the exact function values at the interpolant points. If an empty list is given, i.e., [], depending on the algorithm, a default list of interpolant points is used. Also, the size of the array of ti depends on the algorithm as follows. If kind is:

ext or eig, the no ti is required.
mbf, then a single point ti may be specified.
imbf, spl, or rbf, then a list of ti with arbitrary size may be specified.
crf or rpf, then a list of ti points with even size may be specified.

kwargs**kwargs

Additional arguments to pass to each specific algorithm. See documentation for each kind in the above.

See also

imate.InterpolateSchatten
imate.InterpolateTrace
imate.logdet

Notes

Interpolation of Affine Matrix Function:

This class interpolates the log-determinant of the one-parameter matrix function:

\[t \mapsto \left( \mathbf{A} + t \mathbf{B} \right)^p,\]

where the matrices \(\mathbf{A}\) and \(\mathbf{B}\) are symmetric and positive semi-definite and \(t \in [t_{\inf}, \infty)\) is a real parameter where \(t_{\inf}\) is the minimum \(t\) such that \(\mathbf{A} + t_{\inf} \mathbf{B}\) remains positive-definite.

The interpolator is initialized by providing \(q\) interpolant points \(t_i\), \(i = 1, \dots, q\), which are often given in a logarithmically spaced interval \(t_i \in [t_1, t_p]\). The interpolator can interpolate the above function at arbitrary inquiry points \(t \in [t_1, t_p]\) using various methods.

References

[1]

Ameli, S., and Shadden. S. C. (2022). Interpolating Log-Determinant and Trace of the Powers of Matrix \(\mathbf{A} + t \mathbf{B}\). Statistics and Computing 32, 108. https://doi.org/10.1007/s11222-022-10173-4.

Examples

Basic Usage:

Interpolate the log-determinant of the affine matrix function \((\mathbf{A} + t \mathbf{B})\) using imbf algorithm and the interpolating points \(t_i = [10^{-2}, 10^{-1}, 1, 10]\).

>>> # Generate two sample matrices (symmetric and positive-definite)
>>> from imate.sample_matrices import correlation_matrix
>>> A = correlation_matrix(size=20, scale=1e-1)
>>> B = correlation_matrix(size=20, scale=2e-2)

>>> # Initialize interpolator object
>>> from imate import InterpolateLogdet
>>> ti = [1e-2, 1e-1, 1, 1e1]
>>> f = InterpolateLogdet(A, B, kind='imbf', ti=ti)

>>> # Interpolate at an inquiry point t = 0.4
>>> t = 4e-1
>>> f(t)
2.514109025242562

Alternatively, call imate.InterpolateLogdet.interpolate() to interpolate at points t:

>>> # This is the same as f(t)
>>> f.interpolate(t)
2.514109025242562

To evaluate the exact value of the log-determinant at point t without interpolation, call imate.InterpolateLogdet.eval() function:

>>> # This evaluates the function value at t exactly (no interpolation)
>>> f.eval(t)
2.498580941943615

It can be seen that the relative error of interpolation compared to the exact solution in the above is \(0.62 \%\) using only four interpolation points \(t_i\), which is a remarkable result.

Warning

Calling imate.InterpolateLogdet.eval() may take a longer time to compute as it computes the function exactly. Particularly, if t is a large array, it may take a very long time to return the exact values.

Arguments Specific to Algorithms:

In the above example, the imbf algorithm is used. See more arguments of this algorithm at imate.InterpolateSchatten(kind=’imbf’). In the next example, we pass basis_func_type specific to this algorithm:

>>> # Passing kwatgs specific to imbf algorithm
>>> f = InterpolateLogdet(A, B, kind='imbf', ti=ti,
...                       basis_func_type='ortho2')
>>> f(t)
2.514109025242562

You may choose other algorithms using kind argument. For instance, the next example uses the Chebyshev rational functions with an additional argument func_type specific to this method. See more details at imate.InterpolateSchatten(kind=’crf’).

>>> # Generate two sample matrices (symmetric and positive-definite)
>>> f = InterpolateLogdet(A, B, kind='crf', ti=ti, func_type=1)
>>> f(t)
2.4961565270080617

Passing Options:

The above examples, the internal computation is passed to imate.logdet() function. You can pass arguments to the latter function using options argument. To do so, create a dictionary with the keys as the name of the argument. For instance, to use imate.logdet(method=’slq’) method with min_num_samples=20 and max_num_samples=100, create the following dictionary:

>>> # Specify arguments as a dictionary
>>> options = {
...     'method': 'slq',
...     'min_num_samples': 20,
...     'max_num_samples': 100
... }

>>> # Pass the options to the interpolator
>>> f = InterpolateLogdet(A, B, options=options, kind='imbf', ti=ti)
>>> f(t)
2.5387905608340637

You may get a different result than the above as the slq method is a randomized method.

Interpolate on Range of Points:

Once the interpolation object f in the above example is instantiated, calling imate.InterpolateLogdet.interpolate() on a list of inquiry points t has almost no computational cost. The next example inquires interpolation on 1000 points t:

Interpolate an array of inquiry points t_array:

>>> # Create an interpolator object again
>>> ti = [1e-2, 1e-1, 1, 1e1]
>>> f = InterpolateLogdet(A, B, ti=ti)

>>> # Interpolate at an array of points
>>> import numpy
>>> t_array = numpy.logspace(-2, 1, 1000)
>>> norm_array = f.interpolate(t_array)

One may plot the above interpolated results as follows:

>>> import matplotlib.pyplot as plt

>>> # Plot settings (optional)
>>> from imate._utilities import set_custom_theme
>>> set_custom_theme(font_scale=1.15)

>>> plt.semilogx(t_array, norm_array, color='black')
>>> plt.xlim([t_array[0], t_array[-1]])
>>> plt.ylim([0, 40])
>>> plt.xlabel('$t$')
>>> plt.ylabel('$\mathrm{logdet} (\mathbf{A} + t \mathbf{B})$')
>>> plt.title('Interpolation of Log-Determinant')
>>> plt.show()

Plotting Interpolation and Compare with Exact Solution:

A more convenient way to plot the interpolation result is to call imate.InterpolateLogdet.plot() function.

>>> f.plot(t_array)

In the above, \(f_p(t) = \mathrm{logdet} (\mathbf{A} + t \mathbf{B})\). If you set normalize to True, it plots the normalized function

\[g_p(t) = \mathrm{logdet}(\mathbf{A} + t \mathbf{B}) - \mathrm{logdet}(\mathbf{B}).\]

To compare with the true values (without interpolation), pass compare=True to the above function.

Warning

By setting compare to True, every point in the array t is evaluated both using interpolation and with the exact method (no interpolation). If the size of t is large, this may take a very long run time.

>>> f.plot(t_array, normalize=True, compare=True)

Attributes:

kindstr: Method of interpolation
verbosebool: Verbosity of the computation process
nint: Size of the matrix
qint: Number of interpolant points

Methods

`__call__`(t)	Interpolate at the input point t.
`eval`(t)	Evaluate the exact value of the function at the input point t without interpolation.
`interpolate`(t)	Interpolate at the input point t.
`get_scale`()	Returns the scale parameter of the interpolator.
`lower_bound`(t)	Bound of the interpolation function.
`plot`(t[, normalize, compare])	Plot the interpolation results.