Matrix rank computation

[Mikhail Katliar]

It would be convenient to have a function to compute rank of a matrix:

blaze::DynamicMatrix A(m, n);
// Initialize ...
blaze::size_t r = rank(A);

The rank can be computed as the number of singular values which exceed a specific tolerance. The default tolerance can be computed as max(size(A))*eps(norm(A)), in accordance with how the MATLAB rank() function is implemented: https://www.mathworks.com/help/matlab/ref/rank.html

A draft implementation of the rank() function could look like the following:

/// @brief Computes rank of a matrix
///
/// @param A the whose rank is to be computed.
/// @return Rank of matrix \a A.
template <typename MT, bool SO>
inline size_t rank(Matrix<MT, SO> const& A)
{
    using ET = ElementType_t<MT>;
    DynamicVector<ET> s;

    // Compute singular values of A
    svd(~A, s);

    // Count singular values greater than the tolerance.
    // The tolerance is computed as max(size(A))*eps(norm(A)), 
    // in accordance with how the MATLAB rank() function works:
    // https://www.mathworks.com/help/matlab/ref/rank.html,
    // taking into account that norm(A) == max(svd(A))
    ET const tol = std::max(rows(A), columns(A)) * std::numeric_limits<ET>::epsilon() * max(s);
    return std::count_if(begin(s), end(s), [tol] (auto val) { return abs(val) > tol; });
}

[Klaus Iglberger]

Hi Mikhail!

Thanks for creating this issue and for continually providing great ideas for possible extensions. This indeed sounds like a reasonable extension that we will consider for next releases. Thanks also for the sample code. The solution that we will provide for Blaze will likely look very similar to your code.

Best regards,

Klaus!

[Klaus Iglberger]

Summary

The feature has been implemented, tested, optimized, and documented as required. It is immediately available via cloning the Blaze repository and will be officially released in Blaze 3.7.

The rank() function

The rank() function computes the rank of a given dense matrix:

blaze::DynamicMatrix<double> A( 5UL, 8UL );
// ... Initialization
rank( A );

The rank is determined as the number of singular values greater than a given tolerance. This tolerance is computed as

tolerance = max(m,n) * max(s) * epsilon,

where m is the number of rows of the dense matrix, n is the number of columns of the dense matrix, max(s) is the maximum singular value of the dense matrix and epsilon is the difference between 1 and the least value greater than 1 that is representable by the floating point type of the singular values.

Note that the rank() function 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 is depending on 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.