Provide support for dense matrix inversion

[Klaus Iglberger]

Description

Many engineering applications require the inversion of dense matrices. For this purpose the LAPACK library provides efficient kernels to first factorize/decompose a given matrix and based on the result compute the matrix inversion. However, the interface of LAPACK is very general, not very intuitive, and thus error prone. Consider for instance the following code to invert a given DynamicMatrix in-place:

extern "C" {
void dgetrf_( int* m, int* n, double* a, int* lda, int* ipiv, int* info );
void dgetri_( int* n, double* a int* lda, int* ipiv, double* work, int* lwork, int* info );
}


blaze::DynamicMatrix<double,columnMajor> A( ... );

int m    ( A.rows() );
int n    ( A.columns() );
int lda  ( A.spacing() );
int lwork( n );
int info ( 0 );

std::vector<int> ipiv( blaze::min( m, n ) );
std::vector<double> work( n );

dgetrf_( &m, &n, A.data(), &lda, &ipiv[0], &info );
dgetri_( &n, A.data(), &lda, &ipiv[0], &work[0], &lwork, &info );

The Blaze library should act as a wrapper around the functionality of LAPACK and thus provide an easy and intuitive interface:

blaze::DynamicMatrix<double,columnMajor> A( ... );

invert( A );   // In-place inversion of matrix A
B = inv( A );  // Storing the result of the inversion in B, leaving A unchanged

Tasks

• design the interface for the matrix inversion
• interface with the LAPACK library to provide an efficient and robust inversion algorithm
• provide a full documentation of the feature, including ...
  • ... an explanation of the algorithm
  • ... limitations of the algorithms
  • ... examples of use
  • ... example computations
• ensure compatibility to all existing matrix types
• add appropriate test cases for the feature

[Harry Hill]

This sounds great, do you plan to expose the lapack functions too? In a library for filtering I found it useful to have these functions, https://github.com/hrhill/linalg/blob/master/include/linalg/lapack/blaze.hpp as well as a higher level version with an easier interface. Also, for symmetric pd matrices you could use cholesky to invert rather than lu

[Klaus Iglberger]

Hi!

Thanks a lot for your input. We plan to expose all LAPACK functions both in form of the original interface and via a higher level interface. For instance, these are the signatures for the (P)LU decomposition functions:

namespace blaze {

void getrf( int m, int n, float* A, int lda, int* ipiv, int* info );
void getrf( int m, int n, double* A, int lda, int* ipiv, int* info );
void getrf( int m, int n, complex<float>* A, int lda, int* ipiv, int* info );
void getrf( int m, int n, complex<double>* A, int lda, int* ipiv, int* info );

template< typename MT, bool SO >
void getrf( DenseMatrix<MT,SO>& A, int* ipiv );

} // namespace blaze

These wrapper functions are already available via cloning the Blaze repository and will be officially released in Blaze 2.6. In case you have something else in mind, we are grateful for additional input.

Best regards,

Klaus!

P.S.: Edited to show the final signatures of the wrapper functions.

[Klaus Iglberger]

Summary

The feature has been implemented, tested, optimized (including manual optimization for matrices of up to 6x6) and documented as required. It is immediately available via cloning the Blaze repository and will be officially released in Blaze 2.6.

Dense Matrix Inversion

The inverse of a square dense matrix can be computed via the inv() function:

blaze::DynamicMatrix<float,blaze::rowMajor> A, B;
// ... Resizing and initialization
B = inv( A );  // Compute the inverse of A

Alternatively, an in-place inversion of a dense matrix can be performed via the invert() function:

blaze::DynamicMatrix<double,blaze::rowMajor> A;
// ... Resizing and initialization
invert( A );  // In-place matrix inversion

Both the inv() and the invert() functions will automatically select the most suited matrix inversion algorithm depending on the size and type of the given matrix. For small matrices of up to 6x6, both functions use manually optimized kernels for maximum performance. For matrices larger than 6x6 the inversion is performed by means of the most suited matrix decomposition method: In case of a general or triangular matrix the LU decomposition is used, for symmetric matrices the LDLT decomposition is applied and for Hermitian matrices the LDLH decomposition is performed. However, via the invert() function it is possible to explicitly specify the matrix inversion algorithm:

using blaze::byLU;
using blaze::byLDLT;
using blaze::byLDLH;
using blaze::byLLH;

// In-place inversion with automatic selection of the inversion algorithm
invert( A );

// In-place inversion of a general matrix by means of an LU decomposition
invert<byLU>( A );

// In-place inversion of a symmetric indefinite matrix by means of a Bunch-Kaufman decomposition
invert<byLDLT>( A );

// In-place inversion of a Hermitian indefinite matrix by means of a Bunch-Kaufman decomposition
invert<byLDLH>( A );

// In-place inversion of a positive definite matrix by means of a Cholesky decomposition
invert<byLLH>( A );

Whereas the inversion by means of an LU decomposition works for every general square matrix, the inversion by LDLT only works for symmetric indefinite matrices, the inversion by LDLH is restricted to Hermitian indefinite matrices and the Cholesky decomposition (LLH) only works for Hermitian positive definite matrices. Please note that it is in the responsibility of the function caller to guarantee that the selected algorithm is suited for the given matrix. In case this precondition is violated the result can be wrong and might not represent the inverse of the given matrix!

For both the inv() and invert() function the matrix inversion fails if ...

• ... the given matrix is not a square matrix;
• ... the given matrix is singular and not invertible.

In all failure cases either a compilation error is created if the failure can be predicted at compile time or a std::invalid_argument exception is thrown.

Note that the matrix inversion can only be used for dense matrices with float, double, complex<float> or complex<double> element type. The attempt to call the function with matrices of any other element type or with a sparse matrix results in a compile time error!

Also note that the function inverts the dense matrix by means of LAPACK kernels. Thus the function can only be used if the fitting LAPACK library is available and linked to the executable. Otherwise a linker error will be created.

It is important to note that it is not possible to use any kind of view on the expression object returned by the inv() function. Also, it is not possible to access individual elements via the function call operator on the expression object:

row( inv( A ), 2UL );  // Compilation error: Views cannot be used on an inv() expression!
inv( A )(1,2);         // Compilation error: It is not possible to access individual elements!

This function does not provide any exception safety guarantee, i.e. in case an exception is thrown the matrix may already have been modified.

[Terry Tucker]

Data annotation for specific domains is a game-changer for businesses looking to streamline their operations. The solutions offered by data annotation providers are comprehensive and tailored to meet the specific needs of each industry. Check out this informative article for more information: https://www.technotification.com/2023/06/data-annotation-for-specific-domains-solutions-and-benefits.html. The benefits of using these services are numerous, including increased accuracy, faster turnaround times, and improved data quality. Overall, my experience with data annotation for specific domains has been nothing short of exceptional, and I highly recommend it to anyone looking to optimize their data management processes.