Diagonal matrix from a vector of diagonal elements

[Mikhail Katliar]

I often need to create a diagonal matrix from a vector of diagonal elements. So far my approach is to create a DiagonalMatrix and set its diagonal elements. According to the documentation, "dense diagonal matrices store all elements, including the non-diagonal elements, and therefore don't provide any kind of memory reduction!"

The diagonal matrix created in this way is then multiplied with other dense matrices. I think having a band/diagonal matrix adaptor for a vector would make it not only more convenient, but also more efficient by avoiding unnecessary multiplications with zeros (at least for large enough matrices).

This is probably related to #73.

[Klaus Iglberger]

Hi Mikhail!

Thanks a lot for proposing this idea. It indeed sounds very similar to issue #73, which proposes the introduction of a band matrix. In the context of this issue we will have to think about whether to implement BandMatrix as an adaptor or as a new type of matrix. If we want to support BLAS and LAPACK functionality as suggested in the issue then this will very likely lead to an implementation that is similar to what you are proposing (i.e. a sparse matrix consisting of one or multiple dense vectors). For that reason we see this proposal as a duplicate of the more general idea in issue #73.

However, perhaps Blaze already provides the featuers you need in a different form. First, as stated in the wiki, DiagonalMatrix is an adaptor for both dense and sparse matrices. You can therefore also create sparse diagonal matrices:

// Definition of a compressed row-major single precision diagonal matrix
blaze::DiagonalMatrix< blaze::CompressedMatrix<float,blaze::rowMajor> > D;

From both a performance and memory consumption point of view, this is the preferred choice for large diagonal matrices.

Second, since Blaze 3.3 you can create band views. This allows you to comfortably assign a vector to the diagonal elements of a matrix:

blaze::DynamicMatrix<double,blaze::rowMajor> A;
blaze::DynamicVector<double> x;
// ... Resizing and initialization

// Creating a reference to the diagonal of matrix A via the diagonal() function
auto diag = diagonal( A );

diag = x;  // Assigning vector x to the diagonal of matrix A

Perhaps a combination of a sparse DiagonalMatrix and a band view will work well for you until band matrices have been implemented. Thanks again for creating this issue,

Best regards,

Klaus!

[Klaus Iglberger]

Duplicate of #73.

[Mikhail Katliar]

Hi Klaus,

I still think that having a function which takes a vector argument and returns a matrix expression representing a diagonal matrix whose diagonal is the given vector would be a useful thing.

I can of course create a temporary variable for the matrix:

// Some function taking a matrix argument
template <typename MT, bool SO>
void f(blaze::Matrix<MT, SO> const& m); 

// Some other function
void g(blaze::DynamicVector<double> const& v)
{
  blaze::DiagonalMatrix< blaze::DynamicMatrix<double> > D;
  diagonal(D) = v;  
  f(D);
}

but this is 3 extra lines of code. It becomes even less nice if I need to pass a matrix argument to a ctor in a class member initialization list:

class X
{
public:
  template <typename MT, bool SO>
  X(blaze::Matrix<MT, SO> const& m); 
};

class Y
{
public:
  Y(blaze::DynamicMatrix<double> const& v)
  : x_(makeDiagonalMatrix(v))
  {}

private:
  X x_;

  static auto makeDiagonalMatrix(blaze::DynamicMatrix<double> const& v)
  {
    blaze::DiagonalMatrix< blaze::DynamicMatrix<double> > D;
    diagonal(D) = v;    
    return D;
  }
};

If the requested feature is implemented, it will be just

void g(blaze::DynamicVector<double> const& v)
{
  f(diag(v));
}

class Y
{
public:
  Y(blaze::DynamicMatrix<double> const& v)
  : x_(diag(v))
  {}

private:
  X x_;
};

[Mikhail Katliar]

Or maybe adding another DiagonalMatrix ctor taking a blaze::Vector argument would do the job?

template <typename VT, bool TF>
explicit inline DiagonalMatrix( Vector<VT, TF> const& diagonal_elements );

[Klaus Iglberger]

Hi Mikhail!

I completely agree that the creation of a temporary diagonal matrix in the form you described is not satisfactory. However, if you only need to create the diagonal matrix for a subsequent matrix multiplication, I might have a solution for your problem. Instead of computing a matrix multiplication between a diagonal matrix and a general matrix, you could also perform a Schur product with an expanded matrix. The following code tries to explain the idea:

using namespace blaze;

DynamicMatrix<int> A;     // The general matrix
DynamicVector<int> d;     // The vector containing the diagonal entries
DynamicMatrix<int> B, C;  // The two target matrices 

DynamicMatrix<int> D;  // The diagonal matrix (only required for illustration purposes)

// ... Resizing and initialization

B = A % expand<3UL>( trans(d) );  // Schur product with expanded row vector; Equivalent to B = A * D;
C = expand<3UL>( d ) % A;         // Schur product with expanded column vector; Equivalent to C = D * A;

Via the expand() function (introduced in Blaze 3.5; see the wiki) it is possible to convert a vector into a matrix. In the exceptional case that you need to perform a matrix multiplication with a diagonal matrix, this kind of matrix fits perfectly for a Schur product. An additional benefit is the better performance.

I hope this helps. If not, then we’ll have to think about a simpler way to create a diagonal matrix from a vector.

Best regards,

Klaus!

[Daniel Baker]

By the way, I’ve had a similar need, and I found that the Schur product with expand(v, n) worked correctly, but it seems error-prone in that the user needs to specify the correct other dimension and ensure that it is the diagonal and not the antidiagonal that is used.