imate.AffineMatrixFunction#

class imate.AffineMatrixFunction(A, B=None)#

Create a one-parameter affine matrix function object from input matrices.

Given two matrices \(\mathbf{A}\) and \(\mathbf{B}\), a one-parameter affine operator \(t \mapsto \mathbf{A}(t)\) is created by

\[\mathbf{A}(t) = \mathbf{A} + t \mathbf{B}.\]

The operator accepts various matrix types with a unified interface, establishes a fully automatic dynamic buffer to allocate, deallocate, and transfer data between CPU and multiple GPU devices on demand, as well as performs basic matrix-vector operations for multiple requested parameters with high performance on both CPU or GPU devices.

An instance of this class can be used as an input matrix to any function in imate that accepts slq method.

Parameters:

Anumpy.ndarray or scipy.sparse: The input matrix A can be dense or sparse (both CSR and CSC formats are supported). Also, the matrix data type can be either 32-bit, 64-bit, or 128-bit. The input matrix can be stored either in row-ordering (C style) or column-ordering (Fortran style).

Note

128-bit data type is not supported on GPU.
Bnumpy.ndarray or scipy.sparse, default=None: The input matrix B should have the same size, data type, matrix format, and row-ordering storage as the matrix A. When B is not given (None), it is assumed that B is the identity matrix.

See also

imate.Matrix

Notes

Where to use this class:

The instances of this class can be used just as a normal matrix in any function in imate that accepts slq method. For instance, when calling imate.logdet() function using method=slq argument, the input matrix A can be an instance of this class.

Why using this class:

This class represents an affine matrix function \(\mathbf{A}(t) = \mathbf{A} + t\mathbf{B}\) where \(t\) is a parameter. An instance of this class is suitable to be used in applications where a matrix function of \(\mathbf{A}(t_i)\) (such as its log-determinant, the trace of its inverse, etc) should be computed for a wide range of parameter inputs \(t_i\), \(i = 1, \dots,q\) all at once.

The parameter \(t\) is not set at the instantiation of this class. Rather, a caller function can request a range of matrices \(\mathbf{A}(t_i)\) for multiple values \(t_i\), \(i=1, \dots, q\). This class can provide these matrices one by one or all at once efficiently depending on the case.

Features:

This class internally creates instances of imate.Matrix class for each of the matrices \(\mathbf{A}\) and \(\mathbf{B}\) (if B is not None). This object has the following features.

It acts as a data container that can handle a variety of matrices and data types with a unified interface.
It also creates a dynamic buffer for the data of matrices on either CPU or multi-GPU devices with fully automatic allocation, deallocation, and data transfer between CPU and GPU devices.
It performs basic linear algebraic operations on the matrix data with high performance on CPU or GPU devices.

Further details of each of the points are described below.

1. Unified interface for matrix types:

The linear operator encompasses the following matrix types all in one interface:

Dense or sparse matrices with CSR or CSC sparse formats.

32-bit, 64-bit, and 128-bit data types.

Row-ordering (C style) and column-ordering (Fortran style) storage.

2. Unified interface for memory buffer between CPU and GPU:

This class creates a dynamic buffer to automatically allocate, deallocate, and transfer matrix data based on demand between CPU and multiple GPU devices. These operations are performed in parallel for each GPU device. Also, deallocation is performed by a smart garbage collector, so the user should not be concerned about cleaning the data on the device and the issue of memory leak. An instance of this class can be efficiently reused between multiple function calls. This class uses a lazy evaluation strategy, meaning that data allocation and transfer are not performed until a caller function requests matrix data from this class.

3. Unified interface for basic algebraic operations:

This class handles the following internal operations:

Matrix vector product \(\boldsymbol{y}(t) = \mathbf{A}(t) \boldsymbol{x}\).
Transposed matrix vector product \(\boldsymbol{y}(t) = \mathbf{A}^{\intercal}(t) \boldsymbol{x}\).
Additive matrix vector product \(\boldsymbol{y}(t) = \boldsymbol{y}_0 + \mathbf{A}(t) \boldsymbol{x}\).
Additive transposed matrix vector product \(\boldsymbol{y}(t) = \boldsymbol{y}_0 + \mathbf{A}^{\intercal}(t) \boldsymbol{x}\).

All the above operations are handled internally when an instance of this class is called by other functions.

Each of the above operations has various internal implementations depending on whether the matrix format is dense, sparse CSR, or sparse CSC, whether the memory storage is C style or Fortran style, whether the data type of 32-bit, 64-bit, or 128-bit, whether the operation is performed on CPU or multi-GPU devices, and whether B is identity matrix or generic matrix (to use efficient matrix operations on identity matrix). This class unifies all these implementations in one interface.

Examples

Create an affine matrix function object:

Create a very large sparse matrix with 64-bit data type:

>>> # Create a random sparse matrix with the size of ten million
>>> from imate import toeplitz
>>> n = 10000000
>>> A = toeplitz(2, 1, n, gram=True, format='csr', dtype='float64')
>>> print(A.dtype)
dtype('float64')

>>> print(type(A))
scipy.sparse.csr.csr_matrix

Create the linear operator \(\mathbf{A}(t) = \mathbf{A} + t \mathbf{I}\) from the matrix A where \(\mathbf{I}\) is the identity matrix:

>>> # Import affine matrix function class
>>> from imate import AffineMatrixFunction

>>> # Create a matrix operator object from matrix A
>>> Aop = AffineMatrixFunction(A)

Operation on CPU:

Pass the above linear operator as input to any of the matrix functions in imate that accepts slq method, such as imate.logdet(). Compute log-determinant of \(\mathbf{A}(t)\) for an array of \(t = [-1, 0, +1]\) all at once:

>>> # Import affine matrix function
>>> from imate import logdet

>>> # A list of parameters t to pass to Aop
>>> t = [-1.0, 0.0, 1.0]

>>> # Compute log-determinant of Aop for all parameters t
>>> logdet(Aop, method='slq', parameters=t)
array([ 68.71411681, 135.88356906, 163.44156683])

Operation on GPU:

Reuse the same object Aop as created in the previous example. However, this time by calling imate.logdet() function with gpu=True argument, the data is automatically transferred from CPU to each of the multi-GPU devices (then, it will be available on both CPU and GPU devices).

>>> # Compute log-determinant of Aop on GPU
>>> logdet(Aop, method='slq', gpu=True)
array([ 68.71411681, 135.88356906, 163.44156683])

Attributes:

Anumpy.ndarray or scipy.sparse, (n, n): Input matrix A from python object.
Bnumpy.ndarray or scipy.sparse, (n, n): Input matrix B from python object.
cpu_Aopobject: Matrix object A on CPU.
gpu_Aopobject: Matrix object A on GPU.
gpubool, default=False: If True, the matrix object is created for GPU devices.
num_gpu_devicesint, default=0: Number of GPU devices to be used. If 0, it uses the maximum number of GPU devices that are available.
initialized_on_cpubool, default=False: Indicates whether the matrix data is allocated in CPU.
initialized_on_gpubool, default=False: Indicates whether the matrix data is allocated in GPU.
data_type_namestr, default=None: The type of matrix data, and can be float32, float64, or float128.
num_parametersint, default=1: The number of parameters of the linear operator. For AffineMatrixFunction class, this parameter is always 1 corresponding to the parameter t.

Methods

`initialize`(gpu, num_gpu_devices)	Initializes the object.
`get_num_rows`()	Returns the number of rows of the matrix.
`get_num_columns`()	Returns the number of columns of the matrix.
`is_sparse`()	Determines whether the matrix is dense or sparse.
`get_nnz`()	Returns the number of non-zero elements of the matrix.
`get_density`()	Returns the density of non-zero elements of the matrix.
`get_data_type_name`()	Returns the data type name, which can be float32, float64, or float128.
`get_num_parameters`()	Returns the number of parameters of the linear operator.
`get_linear_operator`([gpu, num_gpu_devices])	Sets the linear operator object on CPU or GPU and returns its object.